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Abstract. High-order derivatives of analytic functions are expressible as Cauchy 
integrals over circular contours, which can very effectively be approximated, e.g., 
by trapezoidal sums. Whereas analytically each radius r up to the radius of conver- 
gence is equal, numerical stability strongly depends on r. We give a comprehensive 
study of this effect; in particular we show that there is a unique radius that mini- 
mizes the loss of accuracy caused by round-off errors. For large classes of functions, 
though not for all, this radius actually gives about full accuracy; a remarkable fact 
that we explain by the theory of Hardy spaces, by the Wiman-Valiron and Levin- 
Pfluger theory of entire functions, and by the saddle-point method of asymptotic 
analysis. Many examples and non-trivial applications are discussed in detail. 
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1. Introduction 

Real variable formulae for the numerical calculation of high-order derivatives 
severely suffer from the ill-conditioning of real differentiation. Balancing approx- 
imation errors with round-off errors yields an inevitable minimum amount of 
error that blows up as the order of differentiation increases (see, e.g., Miel and 
Mooney 1985, Thm. 2). It is therefore quite tricky using these formulae with 
hardware arithmetic, to obtain any significant digits for derivatives of orders, 
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say, hundred or higher. For functions which extend analytically to the complex 
plane, numerical quadrature applied to Cauchy integrals has on various occasions 
been suggested as a remedy (see Gautschi 1997, p. 152/ 187). To be specific, let us 
consider an analytic function / with the Taylor series 1 

CO 

/(z) = J>*z* (\z\<R) (1.1) 

k=0 

having radius of convergence R > (with R = 00 for entire functions). Cauchy's 
integral formula applied to circular contours yields (n — 0, 1,2, . . ., < t < R) 

_/(<<) (0) 

2m J\ z \= r z" +1 
' ' 2n e- ind f{re ie )d6. (1.2) 



2m n Jo 

Since trapezoidal sums 2 are known to converge geometrically for periodic analytic 
functions (Davis 1959), the latter integral is amenable to the very simple and yet 
effective approximation 3 

-1 m— 1 

a n (r,m) = -L- £ e -27djn/m f ^mj/my ^ 
mr j=0 

This procedure for approximating a n was suggested by Lyness (1967). Later, Lyness 
and Sande (1971) observed that the correspondence 



induced by (1.3) is, in fact, the discrete Fourier transform; accordingly they pub- 
lished an algorithm for calculating a set of normalized Taylor coefficients r n a n 
based on the FFT 

Whereas all radii < r < R are, by Cauchy's Theorem, analytically equal, they 
are not so numerically. On the one hand, the geometric convergence rate of the 
trapezoidal sums improves for smaller r. On the other hand, for r — » there is an 
increasing amount of cancelation in the Cauchy integral which leads to a blow-up 
of relative errors (Lyness 1967, p. 130). Moreover, there is generally also a problem 
of numerical stability for r — >■ R (see §3 of this paper). So, once again there arises 
the question of a proper balance between approximation errors and round-off 
errors: what choice of r is best and what is the minimum error thus obtained? 

There is not much available about this problem in the literature. Lyness and 
Sande (1971) circumnavigate it altogether by just considering the absolute errors 
of the normalized Taylor coefficients r n a n instead of relative errors, leaving the 
choice of r to the user as an application-specific scale factor; on p. 670 they write: 



I Without loss of generality, the point of development is z = 0, which we choose for ease of notation 
throughout this paper. Though such series are often named after Maclaurin, we keep the name Taylor 
series to stress that we really do not use anything specific to 2 = 0. 

2 Recall that, for periodic functions, the trapezoidal sum and the rectangular rule are just the same. 

3 For other quadrature rules see the remarks in §2.3. 
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It is natural to ask why this choice of output [i.e., r"a n ] was made, 
rather than perhaps a set of Taylor coefficients a n or a set of deriva- 
tives _p n )(0). The most immediate reason is that the algorithm 
naturally provides a set of normalized Taylor coefficients to a uni- 
form absolute accuracy. If, for example, one is interested in a set of 
derivatives, the specification of the accuracy requirements becomes 
very much more complicated. However, if one looks ahead to the 
use to which the Taylor coefficients are to be put, one finds in 
many cases that uniform accuracy in normalized Taylor coefficients 
corresponds to the sort of accuracy requirement which is most 
convenient. 

Fornberg (1981b, 1981a) addresses the choice of a suitable radius r by suggesting 
a simple search procedure that tries to make (r n a n )™~} approximately proportional 
to the geometric sequence 0.75". If accomplished, this results, for m = 32, in a loss 
of at most about m \ log 10 (0.75)| = 4.0 digits; 4 see §3.1 below. Further, he applies 
Richardson extrapolation to the last three radii of the search process to enhance 
the convergence rate of the trapezoidal sums. However, the success of both devices 
is limited to functions whose Taylor coefficients approximately follow a geometric 
progression. In fact, Fornberg (1981a, p. 542) identifies some problems: 

Some warning about cases in which full accuracy may not be 
reached. Such cases are 

(1) very low-order polynomials (for example, f(z) = 1 + z); 

(2) functions whose Taylor coefficients contain very large isolated 
terms (for example, f(z) = 10 6 + 1 / (1 — z)); 

(3) certain entire functions (for example, f(z) = e z ); 

(4) functions whose radius of convergence is limited by a branch 
point at which the function remains many times [real] differ- 
entiable (for example, /(z) = (1 +z) 10 log(l +z) expanded 
around z = 0). 

As illustrated by the numerical experiments of Figure 1, an answer to the ques- 
tion of choosing a proper radius r becomes absolutely mandatory for derivatives 
of orders of about n = 100 and higher: outside a narrow region of radii there is a 
complete loss of accuracy. However, rather surprisingly, Figure 1 also shows that 
about full accuracy can be obtained for some functions if we choose the optimal 
radius that minimizes the loss of accuracy. We observe that such an optimal radius 
strongly depends on n (and /). This strong dependence, together with the complete 
loss of accuracy far off the optimal radius, prevents us from using, for larger n, 
just a single radius r to calculate all the leading Taylor coefficients «0/ ■ ■ ■ r a n m one 
go; it thus puts the FFT effectively out of business for the problem at hand. 

The goal of this paper is a deeper mathematical understanding of all these 
effects. In particular, we would like to automate the choice of the parameters m 

4\Ve write "=" to indicate that a number has been correctly rounded to the digits given, "~" to 
denote a rigorous asymptotic equality and "as" to informally assert some approximate agreement. 
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Figure i. Numerical stability of using Cauchy integrals to compute 
/("' (0) : plots of the empirical loss of digits (solid red line), that is, the ratio 
of the relative error divided by the machine precision, and its prediction 
by the condition number K(w,r) (dashed blue line) vs. the radius r. The 
vertical lines (dashed green) of the last three plots visualize a finite radius 
of convergence R < oo. In each plot the results for two different orders of 
differentiation are shown: n = 10 (the less steep curves starting from the 
left) and n = 100 (the steeper curves starting farther to the right). The 
number m of nodes of the trapezoidal sum approximation was chosen 
large enough not to change the picture. The qualitative shape (convexity 
in the double logarithmic scale, coercivity and monotonicity properties) 
of these condition number plots can be completely understood from the 
general results in §4. 
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and r and to predict the possible loss of accuracy. This turns out to be a surprisingly 
rich and multi-faceted topic, with relations to some classical results of complex 
analysis such as Hadamard's three circles theorem (§7) as well as to some more 
advanced topics such as the theory of Hardy spaces (§§4/6), the Wiman-Valiron 
theory of the maximum term of entire functions (§8), the Levin-Pfluger theory 
of the distribution of zeros of entire functions (§10); and with relations to some 
advanced tools of asymptotic analysis and analytic combinatorics such as the 
saddle-point method (§9) and the concept of H-admissibility (§11). 

Outline of the Paper. To guide the reader through the thicket of this paper, we 
summarize its most relevant findings: 

• from the point of approximation theory and convergence rates as m — > 00, 
smaller radii are better than larger ones (§2); there are useful explicit upper 
bounds of the number of nodes m in terms of the desired relative error e, 
the order of differentiation n, and the chosen radius r (Eqs. (2.8) and (2.11)); 

• with respect to absolute errors, the calculation of the normalized Taylor 
coefficients r n a n is numerically stable for any radius r < R (§3.1); 

• with respect to relative errors, the loss of significant digits is modeled by 
log 10 /c(«, r) where x(n,r) denotes the condition number of the Cauchy 
integral (§3.2, see also Figure 1), which is independent of the particular 
quadrature rule chosen for the actual approximation; it can be estimated 
on the fly (algorithm given in Figure 3); 

• log k(m, r) is a convex function of logr (Corollary 4.2) and there exists 
an (essentially unique) optimal radius r*(n) = argmin J .K(n, r) that mini- 
mizes the loss of accuracy caused by round-off errors; these optimal radii 
form an increasing sequence satisfying r* (n) — >■ R as n — > 00 (Theorem 4.6); 

• for finite radius of convergence R < 00, the corresponding optimal con- 
dition number k* (n) blows up if f belongs to the Hardy space H 1 (Theo- 
rem 4.7); on the other hand, k* (n) remains essentially bounded if / does 
not belong to the Hardy space H 1 and is amenable to Darboux's method 
(§§5 and 6), in which case there are useful explicit (asymptotic) formulae 
for r*(n) and K*(n) (Eqs. (6.3) and (6.4)); 

• for entire transcendental functions it is more convenient to analyze a cer- 
tain upper bound R(n,r) of the condition number (§7); this yields a unique 
radius r (n) = argmin f , R{n,r), called the quasi-optimal radius, with a cor- 
responding quasi-optimal condition number K (n) = K{n,r {n)) ^ K*(n); 
the quasi-optimal radii also form an increasing sequence with r*(n) — > R 
as n — » 00 (Theorem 7.3); 

• for entire functions of perfectly regular growth there is a simple asymp- 
totic formula for r (n) in terms of the order and type of such a function 
(Theorem 8.4); 

• fo(tt) is the modulus of certain saddle points of |z~"/(z)| in the complex 
plane (Theorem 9.1); the saddle-point method offers a methodology to 
obtain asymptotic results for k (m) (§9.2); 
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• for entire functions of completely regular growth (satisfying certain condi- 
tions on the zeros), the circular contour of radius r (n) is optimal in the 
sense that it passes the saddle points approximately in the direction of 
steepest descent (§10); this yields the extremely simple asymptotic condi- 
tion number bound limsup„ K (n) SC O where O is the number of maxima 
of the Phragmen-Lindelof indicator function of / (Theorem 10.2); in fact, 
there is an explicit asymptotic formula for K<>(tt) in terms of a finite sum 
(Theorem 10.1) that turns out to yield K<>(n) ~ 1 in many relevant examples; 

• for H-admissible entire functions we have K (n) ~ 1 (Corollary 11.3); 

• for entire functions / with non-negative Taylor coefficients the quasi- 
optimal radius r (n) can be calculated as the solution of the scalar convex 
optimization problem r (n) = argmin r r~ n f(r) (Theorem 12.1); we prove 
Ko(tt) ~ 1 for a model of a Fredholm determinant with non-negative Taylor 
coefficients (Eq. (12.8)). 

We shall comprehensively discuss many concrete examples and applications 
throughout this paper: most notably the functions illustrated in Figure 1, the 
functions from the list of the Fornberg quote on p. 3, the functions whose properties 
are listed in Table 2, the functions f(z) = (1 — z)^ (/5 6 R\ No) (Example 5.2), 
the generalized hypergeometric functions (Example 8.2), the reciprocal Gamma 
function f(z) = 1 /T(z) (§10.4), a generating function from the theory of random 
matrices (Examples 3.1 and 12.3), and a generating function from the theory of 
random permutations (Example 12.5). 



2. Approximation Theory 

2.1. Convergence Rates. In this section we recall some basic facts about the con- 
vergence of the trapezoidal sums applied to Cauchy integrals on circular contours. 
We use the notation 

D r = {z G C : \z\ < r}, C r = {z G C : |z| = r}, 

for (open) disks and circles of radius r. Let / be an analytic function as in §1, V m 
be the set of all polynomials of degree ^ m and let 

E m (J;r)= inf ||/-p|| _ (0 < r < R) 

denote the error of best polynomial approximation of / on the closed disk D r . 
Equivalently, by the maximum modulus principle, we have 

E m (f;r)= inf \\f - p\\ L *> { c r ) (0 < r < R). 

The following theorem belongs certainly to the "folklore" of numerical analysis; 
pinning it down, however, in the literature in exactly the form that we need turned 
out to be difficult. For accounts of the general techniques used in the proof see, for 
the aliasing relation, Henrici (1986, §§13.2/4) and, for the use of best approximation 
in estimating quadrature errors, Davis and Rabinowitz (1984, §4.8). 
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Theorem 2.1. Let f be analytic in Dr and < r < R. Then, with the n-th Taylor 
coefficient a n and its approximation a n (r,m) as in (1.2) and (1.3), we have the aliasing 
relation 

r n a n (r,m) = r"'a n i(r,m) (n = n' mod m) (2.1) 
and the error estimate 

r n \a n -a n (r,m)\^2E m _ 1 (f;r) (0^n<m). (2.2) 

Proof. The key to this theorem is the observation that a n (r,m), with ^ n < m, 
is the exact Taylor coefficient of the polynomial p* £ V m -\ that interpolates / in 
the nodes re lm ^ m (j = 0, . . . , m — 1). This fact, and also the aliasing relation, easily 
follows from the discrete orthogonality 

1 n = n' mod m; 

otherwise. 

Now, by introducing the averaging operators 

1 r 2n 1 m-\ 

hif'.r) = f- / e- in6 f{re w )de, Q n {f;r,m) = - £ e ^i«^ f{re lm '' m ), 

A7Z JQ Til j—Q 

(2.3) 

we have r n a n = I n (f;r) and r n a n {r,m) = Q n (/;r,m). The observation about the 
approximation being exact for polynomials implies, for p £ V m -\ and ^ n < m, 
that I n {p',r) = Q n {p;y,m) and hence 

\In(f>r) - Q n (f;r,m)\ 

^ \I„(f;r) - I n (p;r)\ + \Q n (p;r,m) - Q„(f;r,m)\ < 2\\f - p\\i~ {Cr y 
Taking the infimum over all p finally implies (2.2). □ 

From the aliasing relation we immediately infer an important basic criterion for 
the choice of the parameter m, namely the 

Sampling Condition: m > n. (2.4) 

For otherwise, if m ^ n, the value a n (r, m) is just a good approximation of r k ~ n a^, 
with < m the remainder of dividing n by m. However, in general, r^ n a^ will 
differ considerably from a n . 

2.2. Estimates of the Number of Nodes. To obtain more quantitative bounds of 
the approximation error as m — >■ 00, we have a closer look at the error of best 
approximation. With R the radius of convergence of the Taylor series (1.1) of /, 
the asymptotic geometric rate of convergence of this error is given by (Walsh 1965, 
§4-7) 

limsupE m (/;r) 1/m = ^. (2.5) 
Thus, if we introduce the relative error (assuming a n 7^ 0) 



m— 1 



-Inijn I 'm Ircijn' / m 



111 



7=0 



8 



FOLKMAR BORNEMANN 



we get from (2.2) and (2.5) that 



]hnswpS m (n,r) t/m < T -. (2.7) 

m->oo r R 

2.2.1. Finite Radius of Convergence. If R < 00, we obtain from (2.7) that, for n and r 
fixed, 

— \og8 m {n,r)~ l ^ log(£/r) +o(l) (m -> 00). 

Therefore, if m e denotes the smallest value such that S m (n, r) ^ e for m ^ m £ 
(which implies 8 me ~ e as e -> 0), we get the asymptotic bound 

Example 2.2. To illustrate the sharpness of this bound, we consider the function 
f(z) = z/(e z — 1) for n = 100, taking the radius r = 6.22 that is about the optimal 
one shown in Fig. i.e. Here R = 2n and, for a relative error e = 10~ 12 (which is, 
for this particular choice of r, large enough to exclude any finite precision effects 
of the hardware arithmetic), we get 

m £ = 2734 ^ 8 ; p ; -1.00007; 
log(R/r) 

= 2733.80 

thus, the bound (2.8) is an excellent prediction. In Example 6.2 we will see that, 
for general n, the radius r n = 2n(l — n^ 1 ) is, in terms of numerical stability, 
about optimal and yields the estimate m £ w nloge -1 . That is, for e fixed, we get 
m £ = 0(n) as n — > 00, which is the best we could expect in view of the sampling 
condition (2.4). Further examples of this kind are in §§5 and 6. 

2.2.2. Entire Functions. If / is entire, that is, R = 00, the estimate (2.7) shows that 
the trapezoidal sums converge even faster than geometric: 

lim S m (n,r) 1/m =0. 

m— >oo 

In fact, if / is a polynomial of degree d, we already know from Theorem 2.1 that 
the trapezoidal sum is exact for m > d, which implies 5 8 m (n,r) = 0. If / is entire 
and transcendental, a more detailed resolution of the behavior of 5 m depends 
on the properties of / at its essential singularity in z = 00. For example, entire 
functions of finite order p > and type T > (for a definition see §8 below) yield 
(Batyrev 1951, Giroux 1980) 

limsupm 1/ ''E m (f;r) 1/m = riepr) 1 ^. (2.9) 

m— >oo 

We thus get 

)m\suvm l, P5 m {n,rf ,m < r{epr) l/ P (2.10) 

m— foo 

and therefore, for n and r fixed, 
1 1 

— log<5 m (n,r) _1 - -log(m/ (epr)) ^ log(l/r) + o(l) (m -> 00). 



^Recall that we have assumed a„ ^ in the definition of <5„„ which restricts us to n ^ d < , 
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minimal m e jologe 1 /W(loge 1 / 'errP) 



10- 12 

10 -ioo 
-m-1000 



694 



126 



32 48.21 



140.30 



706.73 



Solving for m £ , as defined in §2.2.1, yields the asymptotic bound 



7X77r4^IT77^T( 1 + ( 1 )) ( e ~> °)- 
W(log(e _1 )/(eT7P)) 



(2.11) 



Here W(z) denotes the principal branch of the Lambert W-function defined by 
the equation z = W(z)e w ( z \ In Remark 8.5 we will specify this bound, for entire 
functions of perfectly regular growth, using a particular radius that is about 
optimal in the sense of numerical stability 

Example 2.3. To illustrate the sharpness of this bound, we consider f(z) = e z for 
n — 10 taking the radius r = 10, which we read off from Figure l.a to be close to 
optimal. Here, the order and type of the exponential functions are p = X — 1 (see 
Table 2) and we get the results of Table 1 (that were computed using high-precision 
arithmetic in Mathematica). As we can see, (2.11) turns out to be a very useful 
upper bound. 

2.3. Other Quadrature Rules. To approximate the Cauchy integral (1.2), there 
are other quite effective quadrature rules available besides the trapezoidal sums; 
examples are Gauss-Legendre and Clenshaw-Curtis quadrature. From the point 
of complexity theory, however, Gensun and Xuehua (2005) have shown (drawing 
from the pioneering work of Nikolskii in the 1970s) that the trapezoidal sums 
are, for the problem at hand, optimal in the sense of Kolmogorov. 6 Hence, for 
definiteness and simplicity, we stay with trapezoidal sums in this paper. 

It is, however, important to note that the results of this paper apply to other 
families of quadrature rules as well: first, the estimates (2.8) and (2.11) remain valid 
if the quadrature error is bounded by the error of polynomial best approximation 
(as in (2.2), up to some different constant); which is, e.g., the case for Gauss- 
Legendre and Clenshaw-Curtis quadrature (see Trefethen 2008). Second, the 
discussion of numerical stability in the next section applies to quadrature rules 
with positive weights in general. Ln particular, the estimated digit loss (3.7) depends 
just on the condition number of the Cauchy integral itself, an analytic quantity 
independent of the chosen quadrature rule. Then, starting with §4, optimizing that 
condition number is the main objective of this paper. 



As we have seen in §1 and Figure 1 there are stability issues with using (1.3) in 
the realm of finite precision arithmetic. Specifically, small finite precision errors in 

^That is, the m-point trapezoidal sum minimizes, among all m-point quadrature formulas, the 
worst case quadrature error for the Cauchy integral (1.2) over all analytic functions whose modulus is 
bounded by some constant in an open disk containing |z| < r. 
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s s 

a. gap probability E 2 (10;s) of GUE b. absolute error 

Figure 2. Left: the gap probability E2(10;s) of GUE calculated as the 
10-th Taylor coefficient of a Fredholm determinant; right: the absolute 
error of the calculation. The dotted lines (red) show the results for the 
radius r = 1; the solid lines (blue) show the results for the quasi-optimal 
radius r<>, which depends on s (see Example 12.3 and Figure 7). The 
dashed horizontal lines show the level of machine precision. 

the evaluation of the function / can be amplified to large errors in the resulting 
evaluation of the sum (1.3). This error propagation is described by the condition 
number of the Cauchy integral and depends very much on the chosen radius r 
and on the underlying error concept. 

3.1. Absolute Errors. Any perturbation f of the function / within a bound of the 
absolute error, 

11/ -/1 1 1^(0 ^ e, 

induces perturbations a n {r) and a n {r,m) of the Cauchy integral (1.2) and of its 
approximation (1.3) by the trapezoidal sum. Note that even though the value of the 
Cauchy integral does not depend on the specific choice of the radius r (within the 
range < r < R), the perturbed value a n (r) generally does depend on it. Because 
both the integral and the sum are re-scaled mean values of /, we get the simple 
estimates 

\r n a n — r n a n (r)\ < e, \r n a n (r,m) — r n a n (r,m)\ ^ e. (3.1) 

Thus, the normalized Taylor coefficients r n a n are well conditioned with respect to 
absolute errors (with condition number one); a fact that has basically already been 
observed by Lyness and Sande (1971, p. 670). There are indeed applications were 
absolute errors of normalized Taylor coefficients are a reasonable concept to look 
at, which then typically leads to a proper choice of the radius r. We give one such 
example from our work on the numerical evaluation of distributions in random 
matrix theory (Bornemann 2009). 

Example 3.1. The sequence of hermitian random matrices Xn £ C NxN with 
entries 
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formed from i.i.d. families of real standard normal random variables £h and rjij, 
is called the Gaussian Unitary Ensemble (GUE). 7 The GUE is of considerable interest 
since, on the one hand, various statistical properties of the spectrum ct(Xn) enjoy 
explicit analytic formulas. One the other hand, in the large matrix limit N — >■ oo, by 
a kind of "universal" limit law, these properties are often known (or conjectured) 
to hold for other families of random matrices, too. An example of such a property 
concerns the bulk scaling X^ = 7t~ 1 N 1/ ' 2 Xjv, for which the mean spacing of the 
scaled eigenvalues goes, in the large matrix limit, to one. Basic statistical quantities 
then considered are the gap probabilities 8 

E 2 (n;s) = lim P(#(cr(X N ) n [0,s]) = n), 

N->oo 

the probability that, in the large matrix limit, exactly n of the scaled eigenvalues 
are located in the interval [0, s] . (For Wigner hermitian matrices with a subexpo- 
nential decay, Erdos, Ramirez, Schlein, Tao, Vu and Yau (2009) have, just recently, 
established the universality of £2(0; s).) The generating function of the sequence 
£2(0; s), E2(l; s), E2(2,s), ... is given by the Fredholm determinant of Dyson's sine 
kernel K(x,y)) = sinc(7i(x — y)) (see, e.g., Mehta 2004, §6.4), namely, 
00 . 

£E 2 (k;s)z k =det(l-(l-z)K\ L2{0/S) ). 

k=0 

For given values of n and s, the strategy to calculate £ 2 (w;s) is as follows. First, by 
using the method of Bornemann (2010) for the numerical evaluation of Fredholm 
determinants, the function 

/(z)=det (l-(l- z )K| L2(0/S) ) 

can be evaluated for complex arguments of z up to an absolute error of about 
e = 10~ 15 . Second, the Taylor coefficients E 2 (n;s) of / are calculated by means of 
Cauchy integrals. Now, since these Taylor coefficients are probabilities, the number 1 
is the natural scale for the absolute errors, which makes r — 1 the proper choice for 
the radius (Bornemann 2009, §4.3). By (3.1), we expect an absolute error of about 
e — 10~ 15 , which is confirmed by numerical experiments, see Figure 2. However, 
the figure also illustrates that there is a complete loss of information about the 
tails (that is, those very small probabilities which are about the size of the error 
level or smaller). By controlling the radius with respect to relative errors using the 
method exposed in the rest of this paper, we were able to increase the accuracy of 
the tails considerably. The reader should note, however, that in most applications 
of random matrix theory the accurate calculation of the tails would be irrelevant. 
It typically suffices to just identify such small probabilities as being very small; 
thus the concept of absolute error is completely appropriate in this example. 

There are examples, were small absolute errors of the normalized Taylor coef- 
ficients r n a n are not accurate enough. Because of the super-geometric growth of 

7 In Matlab, the sequence of commands 

X = randn(N) + li*randn(N) ; X = (X+X')/2; 
can be used to sample from the N x N GUE. 

Q 

We denote by #S the number of elements in a finite set S. 
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the factorial, examples of such cases are the derivatives f( n ' (0) = n\ a n , for high 
orders n. Accuracy will only survive the scaling by n! if the Taylor coefficients 
themselves already have small relative errors. 

3.2. Relative Errors. We now consider perturbations / of the function / whose 
relative error can be rendered in the form 



f(re i6 )=f(re i6 )(l + e r (9)), 



< e. 



(3-2) 



Such a perturbation induces a perturbation a n (r ) of the Cauchy integral (1.2) which 
satisfies the straightforward bound (Deuflhard and Hohmann 2003, Lemma 9.1) 



(3-3) 



K{n,r) 



- a„(r)\ 


< jc(n, 


r) ■ e 


a n \ 


7^ 0), where 




Jo 


f(re iS ) 


d6 


rln 

/ e ~ 

Jo 


,ne f(re 


e )d6 



> 1 



(34) 



is the condition number of the Cauchy integral. 9 Note that this number measures 
the amount of cancelation within the Cauchy integral: x(n,r) 3> 1 indicates a large 
amount of cancelation, whereas x(n, r) m 1 if there is virtually no cancelation; see 
Figure 4 for an illustration. 

Correspondingly there are perturbations a n {r,m) of the trapezoidal sum ap- 
proximations (1.3) of the Cauchy integrals. They satisfy the same type of bound, 
namely 

\a„(r,m) -a n (r,m)\ 



\a n {r,m)\ 

of its relative error (assuming a„(r,m) 7^ 0), where 



< K m (n,r) ■ e, 



(3-5) 



m— 1 



E /( 



re 



2mj/m\ 



K m (n,r) 



7=0 



m— 1 



y e -2nijn/m jr^ re ,2nij/r, 
;=0 



^ 1 



(3-6) 



is the condition number of the trapezoidal sum (Higham 2002, p. 538). 

If m is chosen large enough such that the trapezoidal sum a n (r,m) is a good 
approximation of the Cauchy integral a n , then we typically also have 



1 in— 1 

E f(re 2ni ' /m ) 

m r~L 



/=0 



1 

In Jo 



2.T 



f{re m ) 



This is because the integrand \f(re ) \ is a smooth periodic function of 6 and the 
trapezoidal sum therefore gives excellent approximations of this integral, too. 10 



9This condition number is completely independent of how the Cauchy integral is actually computed. 
I0 By the Euler-Maclaurin summation formula, the approximation error is of arbitrary algebraic 
order (Deuflhard and Hohmann 2003, Thm. 9.16). 
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function [val , err , kappa, m] = D(f,n,r) 

fac = exp(ganunaln(n+l) -n*log(r) ) ; 

cauchy = @(t) f ac* (exp(-n*t) . *f (r*exp(t) ) ) ; 

m = max(n+l,8); tol = le-15; 

s = cauchy (2i*pi* (1 :m)/m) ; vail = mean(s) ; errl = NaN; 
while m < le6 
m = 2*m; 

s = reshape([s; cauchy (2i*pi* (1 : 2 :m) /m)] , 1 ,m) ; 

val = mean(s) ; kappa = mean(abs (s) ) /abs (val) ; 

errO = abs (val-vall) /abs (val) ; err = (err0/errl)~2*err0; 

if err <= kappa*tol I I "isfinite (kappa) ; break; end 

vail = val; errl = errO; 

end 

Figure 3. Matlab implementation of calculating (0) using the Cauchy 
integral (1.2) with radius r, approximated by trapezoidal sums. It assumes 
/ to be evaluated up to an relative error tol. The number m of nodes 
is determined by a successive doubling procedure until the estimated 
relative error satisfies a bound corresponding to the level of round-off 
error given by (3.3). The error estimate (see Lyness 1967, Eq. (4.12)) is 
based on the assumption of a geometric rate of convergence (2.5) which 
is excellent if R < 00 and an overestimate if R = 00. The initialization of m 
satisfies the sampling condition (2.4). The doubling of nodes is arranged 
in a way that already computed values of / are re-used. 



Moreover, because of positivity, there are no additional stability issues here. That 
said, for reasonably large m, we have 

K m (n,r) ps K(n,r) 

as long as the computation of a n (r, m) is not completely unstable. We use jc(n, r) in 
the theory developed in this paper; but we use K m (n, r) to monitor stability in our 
implementation that is given in Figure 3. In fact, the examples of Figure 1 show 
that K(n, r) gives an excellent prediction of the actual loss of (relative) accuracy in 
the calculation of the Taylor coefficients; it thus models the dominant effect of the 
choice of the radius r (in fact, for any stable and accurate quadrature rule): 

# lost significant digits ~ log 10 k(m, r). (3.7) 

4. Optimizing the Condition Number 

4.1. General Results on the Condition Number. It is convenient to rewrite the 
expression (3.4) that defines the condition number briefly as 

M x (r) 
\a n \r" 

using the mean of order 1 of the modulus of / on the circle C r , 



(4-2) 
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Concerning the properties of M\ we recall the following classical theorem, for the 
standard proof see Dienes (1931, p. 156) or Polya and Szego (1964, III.310). 

Theorem 4.1 (Hardy 1915). Let f be given by a Taylor series with radius of convergence 
R. The mean value function M\ satisfies, for < r < R: 

(a) Mi (r) is continuously dijferentiable; 

(b) iff ^ 0, log Mi (r) is a convex function of log r; 

(c) if f ^ const, M\{r) is strictly^ increasing. 

Because of log x(n, r) = log Mi (r) — log \a n \ —n log r, there are some immediate 
consequences for the condition number. 

Corollary 4.2. Let f ^ be given by a Taylor series with radius of convergence R. Then, 
for n with a n 7^ and for < r < R: 

(a) K{n,r) is continuously differentiable with respect to r; 

(b) log x(n, r) is a convex function of log r. 

We now study the behavior of n{n, r) as r — > and r — > 00. The first direction is 
simple and gives us the expected numerical instability for small radii. 

Theorem 4.3. Let f ^ be given by a Taylor series with radius of convergence R and let 
a„ be its first non-zero coefficient. Then, for n > no, 

x(n, r) — > 00 (r — > 0); 

but k(«0/ t) — > 1. 

Proof. From the expansion 

Mi(r) = ^ \f( re ' 6 )\ de = \a„ \r n ° + O(r"o +1 ) (r -> 0) 

we get 

K(n,r) ~ (r ->■ 0) 

which implies both assertions. □ 

The other direction, r — » R, is more involved and depends on further properties 
of /. Let us begin with entire functions (R = 00). 

Theorem 4.4. Let f be an entire function. If f is transcendental then, for all n G N, 

K,(n,r) — > 00 (r — > 00). 
If / /s a polynomial of degree d then this results holds for all n =/= d, but x(d, r) — > 1. 
Proo/. Let us assume that, for a particular meN, 

liminf = i imin f _L_ f 2n \f(re ie )\de < 00. 

Then, for all n > m, 

0^ \a n \ ^liminf— \f{re l9 )\ d0 = 0; 

II The fact that the monotonicity is strict has been added to Hardy's theorem by Taylor (1950). 
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that is, a n — 0; implying that f is a polynomial of degree d m. This proves the 
assertion for transcendental /; and for the cases n < d if / is a polynomial of 
degree d. The cases n > d follow trivially from a n = which implies K(n,r) = oo. 
Finally the case n = d gives, because of |/(z)| = \a d \ \z\ A + 0{\z\ d - 1 ) as z — > oo, 

K (d,r) = — i — i f 2n \f(re ie )\d6 = l + 0(r- 1 ) (r -> oo), 
2n\a ( j\r a Jo 

which completes the proof. □ 



For finite radius of convergence, R < oo, we recall the definition of the Hardy 
norm (the last equality follows from the monotonicity of Mi): 

II/IIhi(d r ) = SU P M lW = HmM^r). (4.3) 

If ||/||h1(Dk) < 00 the function / belongs to the Hardy space H 1 (D K ). From the 
strict monotonicity and differentiability of M\ (r) we infer that the function 

cr(r) = log Mi (r) 

satisfies cr'(r) > (0 < r < R). Since log Mi (r) is convex in logr, the function 
ro" '(r) is monotonically increasing. Therefore, the limit 

v = sup rcr 1 '(r) = lim ra'{r) > (4.4) 



0<r<R 



exists (with v = 00 a possibility, however). 

Theorem 4.5. Let f be given by a Taylor series with finite radius of convergence R < 00. 
Then, for a n 7^ 0, 

Mm *(",')= "{" Hl f R) - 

T/zz's is finite if and only if f belongs to the Hardy space H 1 (Dr). Ifn>v then x(n,r) 
is strictly decreasing for < r < R; whereas if v = 00 then, for all n, Kin, r) is strictly 
increasing in the vicinity ofr = R. 

Proof. The limit can be directly read-off from (4.3). If n > v, we have 

-^logzc(?z, r) = cr'(r) - nr^ 1 ^ (v - n)r~ l < (0 < r < R), 

which shows that x{n,r) is strictly decreasing. If v = 00 then cr'(r) — >■ 00 as r — >■ K, 
which implies 

— - log K.{n,r) = cr'(r) — nr~ x — >■ 00 (r — )■ R). 
Hence, k(h 7 r) must be, for r close to K, strictly increasing. □ 
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4.2. The Optimal Radius. Optimizing the numerical stability of the Cauchy inte- 
grals means, by (3.7), to choose a radius r that minimizes the condition number 
x(n,r). The general results of §4.1 imply that such a minimum actually exists. In- 
deed, assuming n > Uq (see Theorem 4.3), a n 7^ 0, and that / is not a polynomial, 12 
we have the following ingredients allowing the optimization: 

• continuity: x(n,r) is continuous for < r < R (Corollary 4.2) and, if R < 00 
and ||/||h1(d r ) < °°' can ^ e continuously continued to r = R (Theorem 4.5); 

• convexity: logje(n,r) is convex in logr (Corollary 4.2); 

• coercivity: x(n,r) — >■ 00 as r — >■ (Theorem 4.3) and, if R = 00 (Theorem 4.4) 
or if R < 00 and H/HhVDr) = 00 (Theorem 4.5), as r — >■ R. 

Hence, by the strict monotonicity of the logarithm, the optimal condition number 

K*(n) = min ;c(n,r) (4.5) 
0<i%R 

exists and is taken for the optimal radius 13 

r*(n) = argminjc(n,r). (4.6) 

0<r^R 

Because the functions r~ n M\(r) and %{n,r) just differ by a factor that is indepen- 
dent of r (namely, \a n \), it is convenient to extend the definition of the optimal 
radius r*(n) to the case a n — by setting 1 ^ 

r*(n) = argminr~"Mi (r). (4.7) 

0<r<R 

Theorem 4.6. Lei f/ze non-polynomial analytic function f be given by a Taylor series with 
radius of convergence R. Then, the sequence r*(n) satisfies the monotonicity 

r*(n) < r*(n + 1) (n > no) 

and /las the limit lim„^ ro r t . (n) = R. Furthermore, the case (n) = R is characterized 
by 

r* (n) = R K < 00, ||/|| h1(Dr) < 00, and i/ < 00, 

and 

R < 00, ||/|| h i( Dr ) < 00, and v < n r*(n) = R. 

Proo/ Because of the optimality of r*(n) and since Mj(r) > 0, we have, for < 

r < r*(n), 

r,(n)-(" +1 )Mi(r*(n)) < r.(n)- 1 r- B Mi(r) < r-^+^M^r). 



I2 Polynomials are addressed by Theorem 4.4: First, one detects the degree d from lim MCO K(rf, r) = 1; 
then, the cases n < d are dealt with as for entire transcendental / of order p = (see §8). 

I3 Since we have no proof of strict convexity, we cannot exclude that the minimizing radius happens 
to be not unique (even though we have not encountered a single such example). However, because of 
convexity, the set of all minimizing radii would form a closed interval. We therefore define r t (n) as the 
smallest minimizing radius; which, in view of (2.7) and (2.10), gives the best rates of approximation of 
the trapezoidal sums. 

I 4Note that all the qualitative results that we stated in §4.1 for %(n,r) hold verbatim for r~ n M\(r), 
independently of whether a„ 7^ or not. 
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Hence, the optimal radius r*(n + 1) must satisfy r*(n + 1) ^ r*(n). This mono- 
tonicity implies that ;'o = lim„^oo r* [n) exists. Let us assume that r$ < R. Then, 
for each r$ < r < R, by taking the limit n — > oo in 

u{n)- l M l {n{n)) x/n ^ r- 1 M 1 (r) 1/n , 

and recalling the continuity of Mi, we conclude r^ 1 ^ r _1 . Since this contradicts 
the choice tq < r, we must have 1'q = R. The characterization of r* (n) = R follows 
straightforwardly from Theorem 4.5. □ 

Bounded analytic functions / that belong to the Hardy space H 1 (D#) are known 
to possess boundary values (Garnett 1981, §11.3); that is, the radial limits 

f(Re ie ) = lim f(re w ) 
exist for almost all angles 8. These boundary values form an L 1 -function, 

\\f\\HHD R) = ^J 2n \f(^)\de r 

whose Fourier coefficients are just the normalized Taylor coefficients of /: 

1 rln 

a n R n = ^ e- me f(Re ie )d6 (n = 0,1,2,...). 



27T Jo 

As the following theorem shows, this fact is bad news for the optimal condition 
number of such functions for large n: it grows beyond all bounds, at a rate that is 
all the more faster the more regular the boundary values of / are. 

Theorem 4.7. Let the analytic function f be given by a Taylor series with finite radius of 
convergence R < 00. If f g H^Dr) then 

k* (n) — > 00 (n — > 00). 

For boundary values of f belonging to the class 15 C kA (C^) the optimal condition number 
grows at least as fast as K*(n) ^ cn k+a for some constant c > 0. 

Proof. Since a n R n are the Fourier coefficients of the L 1 -function formed by the 
radial boundary values of /, the Riemann-Lebesgue Lemma implies 

a n R n — >■ (n — > 00); 

with a rate 0(n~ k ~ a ) if these boundary values belong to the class C k,a (see, e.g., 
Zygmund 1968, §11.4). By Theorem 4.6 we have r*(n) — > R. Hence, for n — >■ 00, 

*.(») = K(n,r,(n)) = fiM")) „ ^ U/lltf(D») , „ 



|«nk*(")" |«nk*(«) n |«n|R" 
since |/|| h j (d r ) > ^ (otherwise we would have / = and _R = 00). □ 



I 5c'c,« denotes the functions that are k times continuously differentiable with a fc-derivative satisfying 
a Holder condition of order < a < 1. 
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5. Examples of Optimal Radii 

Qualitatively, the general results of Section 4.1 are nicely illustrated by the 
examples of Figure 1. In this section we study a couple of important examples 
more quantitatively for large n. 

Example 5.1. This example illustrates the excellent behavior of certain entire 
functions; a general theory will be developed in §§7-12. Here, we consider one of 
the simplest such functions, namely the exponential function 

/CO = 

which is an entire function (R = 00) with the Taylor coefficients a n = l/n\. The 
mean value of the modulus is explicitly given in terms of the modified Bessel 
function of the first kind of order zero (Watson 1944, §3.71), 

Mi(r) = ^ fj I exp(re ie )\d8 = ± j\ rcose d6 = I (r) (r > 0). 

Hence, the condition number is 

K(n,r) = r~ n n\ Io(r). 

Figure 4 illustrates the vast cancelations that occur in the Cauchy integral for large 
condition numbers K.(n, r), that is, for far-from-optimal radii r. Using Stirling's 
formula and the asymptotic expansion of the modified Bessel function (Andrews, 
Askey and Roy 1999, Eq. (4.12.7)), 

we get an explicit description of the optimal radius and its condition number: 
namely, as n — » 00, 

r*{n) = n+ 1 - + ^ i +0{n- 2 ), (5.1a) 

K * (H) = 1 + ri + 28^ +0 ^ 3) - (5 ' lb) 

In fact, already the first term of this expansion for r*(n) gives uniformly excellent 
condition numbers: 

l<x(n,n) < 1.3 (n>l). 

Thus the derivatives of the exponential function can be calculated to full accuracy 
using Cauchy integrals, for all orders n. On the other hand, Figure l.a shows that, 
by choosing a fixed radius r independently of n, it would be impossible to get 
condition numbers that remain moderately bounded for orders of differentiation 
between, say, 1 and 100. This explains the failure that Fornberg (1981a, p. 542) has 
documented using his implementation for the exponential function. 

Example 5.2. In preparation of §6 we consider the family 

/p(z) = (l-z)" ()6eR\No) 
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Re(,-7K*)) &(,-■/('<")) &(,--/(,,")) 




-3-2-1 12 3" 

a. r = 1, K(n,r) = 1.182 x 10 158 b. r = 100, K(n,r) = 1.002 c. r = 200, K(n,r) = 1.502 x 10 : 



Figure 4. Real part (oscillatory, blue line) and absolute modulus (enve- 
lope, red line) of the integrand of the Cauchy integral (1.2) for various 
radii r; f(z) = e z , n = 100. Clearly visible is the huge amount of cancela- 
tion if the condition number k(u, r) is large. Note that this is not an issue 
of frequency, which is moderate and perfectly dealt with by the sampling 
condition (2.4), but rather an issue of amplitude. 



of analytic functions, which are not polynomials for the values of /5 considered. 
The radius of convergence of the Taylor series is R = 1 and the Taylor coefficients 
are given by 



n-f>-\ 



(n = 0,1,2,...). 



By a simple transformation of Euler's integral representation (Andrews et al. 1999, 
Thm. 2.2.1), the mean value of the modulus can explicitly be expressed in terms of 
the hypergeometric function 2F1 (a,b;c;z): 



Mi(r) 



2n Jo 



2/T 



\l-re w fdt 
(1+ryVi 



In 

2/r Jo 
1 



\/l +r 2 -IrcosQ 
4r 



2' 2 ; 1; (1 + r) 2 



(0 ^ r < 1). (5.2) 



The classical results of Gauss (Andrews et al. 1999, Thms. 2.1.3/2.2.2) about the 
hypergeometric function 2F\(a,b;c;z) as z — > 1 imply, as r — >■ 1 from below, 



M 1 {r) 



0rr(g + i) 



jS+1 
2 



2^77 r 



(1-r) 



3+1 



(]6>-i); 

(jS = -1); 
(]6<-l). 



(5-3) 



Therefore, we have to distinguish three cases. 
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Case I: f> > — 1. Here, (5.3) implies that /g belongs to the Hardy space H 1 (Di) 
with norm 

WhHd^ = UmMiCr) = - ) \ > 1. 

^ v^r(J + ij 

(The estimate from below follows from the fact that ||/^||h1(Di) ^ s a convex an d 
coercive function of /3, taking its minimum at p = 0.) The constant v, defined in 
(4.4), can be computed from 



M[(r) = 0(1 + r)"- 3 ((1 +r) 2 2 F X 1; ^ [ |-^) 



Kr-i) 2 fi('ii-?;2; 4r 



2' 2' ' (1 + r) 2 



to have the value 



lim rcr(r) - 



r-H w Mi(l) 2' 
Thus, by Theorems 4.5 and 4.6, the condition number K(n,r) is strictly decreasing 
for n > /5/2 (see Figure l.f for an example); hence 

r*(n) = l (n>0/2), (5.4a) 

which induces (by Stirling's formula) 



**(n) 



^ IIh1(Di) > 1 |r(-p)|n0 +1 -x» (n ->«>). (5.4b) 



m— j8— 1 



This means that for each radius r there will be a complete loss of digits for n 
large enough (e.g., there is already a more than 12 digits loss for f> = 11/2 and 
n = 100, see Figure l.f); an effect that will be the more pronounced the larger /5 is. 
Note that a larger /3 corresponds to higher order real differentiability at the branch 
point 2 = 0; an observation which is in accordance with Theorem 4.7 and which 
helps to explain the failure that Fornberg (1981a, p. 542) has documented using 
his implementation for such functions. 

Case II: f> = — 1. Now, (5.3) shows that /« does not belong to the Hardy space 
H 1 (Di) anymore. Thus, by Theorems 4.5 and 4.6, we have < r*(n) < 1 with 
r*(n) — >■ 1 as ft — >■ 00. Because of a„ = 1, and by (5.3) once more, there is the 
asymptotic expansion 

. , Mi(r) 1 , / 1 \ 
K(n,r) = iW ~ log (r^l). 

It is now a more or less straightforward exercise in asymptotic analysis (de Bruijn 
1981, Chap. 2) to get from here to the following expansions of the optimal radius 
and condition number: as n — » 00, 



r*{n) 



1 ~ / log log n\ , 
nlogn \ft(logn) z / 

K ,(„) = lo | M +0(loglogn). (5.5b) 
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This logarithmic growth is very moderate; indeed, one has 

l< K (n,l - 1 — )<4.8 (3^n^ 10000), 

\ n log n J 

which means that less than one digit is lost for a significant range of n. 

Remark 5.3. In practice it is not always advisable to use the optimal radius: a 
small sacrifice in accuracy might considerably speed up the approximation of the 
Cauchy integral by the trapezoidal sum. In fact, if we recall (2.8), we realize that 
the near-optimal choice r n = 1 — (nlogn) -1 would need about 16 

m £ « n log n ■ log e _1 (5.6) 

nodes to achieve an approximation of relative error e. We can actually get rid 
of the factor logn here if we use the sub-optimal radius f„ = 1 — cm -1 (a > 0) 
instead. Asymptotically, as n — >■ 00, the condition number is then 

K(n,? n ) ~ ^ log = „ {1 _\ n - 1)n tos(«/*) ~ "kg* (5-7) 

and therefore still of logarithmic growth: compared to r n we additionally sacrifice 
just about log 10 e u = 0.43 a digits, independently of n. However, the corresponding 
number of nodes now grows like 

m £ k, ^loge" 1 , (5.8) 

which is about an oc log n improvement in speed. 

To be specific, let us run some numbers for n = 100: Since k(100, fioo) = 
3.25, we are about to lose 0.51 digits using r n ; in hardware arithmetic we could 
therefore strive for a relative error of e = 2 x 10~ 15 . By (5.6) we have to take 
about m £ ps 16000 nodes; actually, a computation with m = 20000 gives us the 
relative error 2.6 x 10~ 15 . In contrast, for a. = 4, we have k(100, f\oo) = 101.63, so 
we are about to lose 2.0 digits using ?„; we could therefore strive for a relative 
error of e = 5 x 10 -14 here. Because of (5.8) we now have to take just about 
m £ ps 800 nodes; and indeed, a computation with m = 800 gives us the relative 
error 4.9 x 10 -14 . Thus, sacrificing just a little more than one digit cuts the number 
of nodes by a factor of 25 (the prediction was 4 log 100 = 18.4). 

Case III: f> < —1. As for f> = —1, (5.3) shows that these fp do not belong to the 
Hardy space H 1 (Di). Thus, by Theorems 4.5 and 4.6, we have < r* (n) < 1 with 
r* (n) — >■ 1 as n — > 00; hence, (5.3) implies the asymptotic expansions 

6 + 1 

r*(n) = l + t- hO(n" 2 ) (n ->■ 00) (5.9a) 



Note that, by (2.8) and (2.11), estimates of the form m e <s • • • include, among other approximations, 
a factor of the form 1 + 0(1) as e — » 0. Therefore, one should not expect too much precision of such 
estimates, in particular not if additionally finite precision effects come into play for e close to machine 
precision. Even then, however, in all the examples of this paper, we observe ratios of the actual values 
of m e to their estimates that are smaller than 1.3; thus, these rough estimates are, in practice, quite 
useful devices to predict the actual computational effort. 
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of the optimal radius and 

i r(-m(-^Y +1 

K*(n) 



2^ 



2 



11 



c^-^-iiVi^' (5 . 9b) 



7T 

of the optimal condition number. Note that there is no explosion in n and that 
Cp — S> 1 monotonically from above as /3 — >■ — oo. Quantitatively we have 

1 ^ c p ^ 2 (/J < -1.362), 

that is, we are just about to lose one binary digit of accuracy within this range of 
values of /3 (for large n). Finally, to accomplish an approximation of relative error 
e by using a trapezoidal sum, we would need, in view of (2.8), about the following 
number of nodes: 

m e w _o_ 1 ■ loge" 1 . (5-10) 

Here are some actual numbers: for f> = —6,n = 100, r„ = l + (/3 + l)n~ 1 , and 
the accuracy requirement e = 10~ 15 , we get 

K(n,r n ) = 1.0769, m e w 700. 

In fact, a computation in hardware arithmetic secures a relative error of 4 x 10 -15 
using m = 900 nodes. 

Example 5.4. We analyze a further example that Fornberg (1981a, p. 542) has 
documented to fail his implementation: 

/(z) = (l+z) 10 log(l+z) 

with radius of convergence R = 1. Having norm H/HjfUrjj) = 180.14, this function 
belongs to the Hardy space H 1 (Di). Theorem 4.7 gives K*(n) — »■ 00 as n — > 00. 
More quantitatively we get, by Theorem 4.5, 

180 14 

K(n,r) ^ k(w,1) = , '-. (n >v = 5.727). 
|«n| 

The asymptotics (the first equality is valid for n ^ 11) 
_ (-1)"" 1 (-l) n_1 10! 



" 11(A) 

implies 



(n — > 00) 



K(n, r) ^ K(n, 1) = 1981.57 1 *\ ~ 5.46 x 10 -4 ■ n 11 (n ->■ 00). 

For instance, n = 50 gives k(50, r) ^ k(50, 1) = 7.4 x 10 13 ; meaning that a loss of 
more than about 14 digits is unavoidable here. 
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Example 5.5. The final example of this section is also taken from the list of failures 
documented by Fornberg (1981a, p. 542): 

/(Z ) =10 6 + _L_ 

with radius of convergence R = 1. This function is a perturbation of the function 
f-i from Example 5.2. Denoting by M\{f^\;r) the mean value of the modulus of 
/_l we get, using (5.3), 



Mi(r) ^ 10 6 + Mi(/li;r) ~ 10 6 + — ^ ' V (r 



71 

,-1 



The sub-optimal choice r n = 1 —n (see Remark 5.3) yields 
K(n,r n ) ^ 10 6 e+ -logn » 3 x 10 6 (1 < n < 10 



7T 



Hence, we expect a loss of (at most) about 6.5 digits throughout this huge range 
of n. The estimate is, in fact, quite sharp: for instance, n = 100 yields 

K(W0,r m ) =2.7 x 10 6 . 

An actual calculation using a trapezoidal sum with m = 4096 nodes yields a 
relative error of 3.13 x 10 -10 which corresponds to a loss of a little more than 6 
digits in hardware arithmetic. 

6. Functions Amenable to Darboux's Theorem 

Example 5.2 contains, in fact, all the information that is needed to address a 
large class of analytic functions: 

/(z) = (l-z)Vz) (/$eR\N ), 

where o(z) is analytic in a neighborhood of Di, z?(l) 7^ 0. In particular, the radius 
of convergence is R = 1. By Darboux's theorem (Wilf 2006, Thm. 5.3.1), the Taylor 
coefficients are asymptotically given by 

a » ^( 1 )f^y( 1 +°( n ~ 1 )) (n->°°). (6.1) 
Hence, the condition number is asymptotically described by 

The mean value of the modulus satisfies, as r — > 1, (compare with (5.3)) 

(j6 > -i); 



Mi(r) = i \l-re i0 f ■ \v{re ie )\d9 ~ < 



27r7o 



[c(l-r)^ 1 (]6<-l). 
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Here, c denotes some positive constant that depends on v and /3. This implies, 
such as in Example 5.2, that, asn^oo, 



= 1 (6 > -1); 



r*(n) < 




(j6 > -1); 



l ~^o%n ( ^ = " 1); *.(»)~ <jclogn (j6 = -l); (6.3) 
6 + 1 U (6<-l). 

n 

For large orders of differentiation, this means that, once more in accordance with 
Theorem 4.7, the Hardy space case 6 > — 1 yields polynomial growth of the 
condition numbers; whereas for 6 = — 1 we get just logarithmic growth and for 
6 < — 1 there is a uniform bound of the condition number. 

To address the last two cases more quantitatively, we can estimate the mean 
modulus by 

Mi(r)< ^^°°( p i) [ 2n \i-re l9 fd9 (0 < r < 1) 

Z7T JO 

with the help of yet another Hardy space norm, defined by 

II/IIh~(D,.) =esssup|/(re' e )|. 

Denoting the condition number of the Cauchy integral for the function fp by 
K(fp;n,r) (recall that this expression can be evaluated in terms of the hypergeo- 
metric function, see (5.2)), we thus obtain a useful estimate of the condition number 
itself, namely 

K{ - H ' r) ^ M \v{l)\ X) <fc n ' r ^ (° <'<!)• ( 6 -4) 
Note that there is nothing special about R = 1 here. For functions of the form 

f{z) = (zQ-z)h{z) (/3eR\N ) 

with |zo| = R, v(z) analytic in a neighborhood of Dr, and v(zq) 7^ we get 
accordingly 

<n,r) ^ ^fe^ K(fp-n,r/R) (0 < r < R). (6.5) 

If there is more than one singularity on the circle Cr, we would have to use sym- 
metry arguments or we would have to consider superpositions of these estimates. 

Example 6.x. We study the example of Figure l.d, that is, 

f(z) = sec(z) 6 

which has radius of convergence R = 71/ 2. To begin with, we extract the poles at 
z = ±7i/2 by the factorization 

f(z) =g{zf-v{z), v{±n/2) = l, 

with the rational function 

. , res„ /o sec res „ /t sec 4tt 2 
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One easily checks that ||p||h°°(Di) = 1/ so that, by (6.5) and by a symmetry argu- 
ment, 

1 < K (n,r) <: K {f„ 6 ;n,r/R) (R = zr/2). 
In view of (6.3) we choose the radius 

n (_ 5" 



" 2 

and obtain (see (5.9b) for a definition of eg) 

_i 9e 5 

1 < K(n,r„) < jc(/_ 6 ;n,l - 5n ~ c_ 6 =——= 1.0686 (n -> 00). 

IzdU 

We should thus be able to get about full accuracy for large orders of differentiation. 
In fact, for n = 100, we have 

K(n,r n ) = 1.0767 < 1.0769 = jc(/_ 6 ; n,r„). 

Striving for a relative error of e = 10 -15 requires, see (5.10), a trapezoidal sum 
with a number of nodes of about 

m £ « I log e" 1 « 700. 

In fact, an actual computation with m = 880 yields a little more than 14 correct 
digits in hardware arithmetic. 

Example 6.2. In this example we address the accurate computation of the Bernoulli 
numbers Bj- given by their exponentially generating function (see Figure i.e) 

which has radius of convergence R = In. We extract the poles at z = ±27B by the 
factorization 

f(z)=g(z)>v(z), v(±2m) = l, 

with the rational function 

= res 27r ;/ res_ 2 7Ti/ = 8tt 2 
*^ j z-2m z + 2m 47r 2 + z 2 ' 

One easily checks that |MIh»(Di) = 27r/(l - e^ 171 ) = 6.2949, so that, by (6.5) and 
by a symmetry argument, 

1 < je(n, r) < 6.3 k(/_i; n, r/i?) (2? = In). 

Because of (5.7) we expect just a moderate loss of accuracy using the choice 
r„=27r(l — n^ 1 ). In fact, for n = 100 we get k(100, rioo) = 7.2355, meaning a loss 
of less than one digit. In view of (5.8) we expect to accomplish an approximation 
error e — 10~ 15 using a trapezoidal sum with a number of nodes of about 



m e 



nloge" 1 » 3500. 



In fact, an actual calculation with m = 4096 gives more than 15 correct digits in 
hardware arithmetic. 
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7. The Quasi-Optimal Radius 

For entire transcendental functions, it turns out that an upper bound of the 
condition number is actually easier to analyze, namely 

*(n,r) = ^<^ = *(»,r), 
a„ \r n \a„\r n 



where 



Mir) = max \f(re' e )\ 



denotes the maximum modulus function of /. In fact, we will see in §§9-12 that 
the radius that is optimal for this upper bound is in many cases already close to 
optimal for the condition number itself. 

For the maximum modulus, the analogue of Hardy's theorem is a classical 
theorem of complex analysis (the three circles theorem); for the standard proof see 
(Markushevich 1977, Vol. II, p. 221) or Polya and Szego (1964, III. 304/305): 

Theorem 7.1 (Hadamard 1896, Blumenthal 1907, Faber 1907). Let f be given by a 
Taylor series with radius of convergence R. The maximum modulus function M satisfies, 
forO<r< R: 

(a) M(r) is continuously differentiate, except for a set of isolated r; 

(b) iff(z) is not a monomial, logM(r) is a strictly convex function of log r; 

(c) if f ^ const, M(r) is strictly increasing. 

With the same proofs as in §4.1 for the condition number, we deduce from 
this theorem the following results (restricting ourselves to entire transcendental 
functions, though). 

Theorem 7.2. Let f be an entire transcendental function with Taylor coefficients a n , and 
let a„ Q be its first non-zero coefficient. Then, for n > uq, with a n 7^ and r > 0: 

(a) x(n, r) is continuously differentiable, except for a set of isolated r; 

(b) log R(n, r) is a strictly convex function of log r; 

(c) ic(n,r) — > 00, as r — >• and r — > 00. 

The same reasoning as in §4.2 shows the existence of the optimal upper bound 

iCo(n) = minic(n,r), (7.1) 

which is now taken for the radius 

r (n) = argminic(n,r). (7.2) 

Note that r (n) is unique because of the strict convexity stated in Theorem 7.2. As 
for r* (n) it is convenient to extend the definition of r (n) to the case of a n = by 
setting 

r (n) = argminr~"M(r). (7.3) 

r>0 

We call r (n) the quasi-optimal radius and define, accordingly, the quasi-optimal 
condition number by 

Ko(n) = K(n,r (n)) ^ K*(n). (7.4) 
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Finally, by repeating the proof of Theorem 4.6 we get: 

Theorem 7.3. Let f be an entire transcendental function. Then, the sequence r (n) 
satisfies the monotonicity 

r<>{n) < r (n + 1) (n > n ) 
and has the limit lim,,-^ r (n) = 00. 

It turns out that the radius r (n) is generally much easier to calculate than the 
optimal radius r*(n) (see Theorems 8.4 and 9.1). Surprisingly in all of these cases 
the radius r (n) is also very close to optimal and the condition number k (m) is 
close to one. Before giving a theoretical frame for these effects, we illustrate them 
by two examples. 

Example 7.4. Since its Taylor coefficients are positive, the exponential function 
f(z) = e z has the maximum modulus function M(r) = e r . A short calculation 
shows that 

r (n)=n, ic<>(n) = ft! = \ / 2nn (1 + 0(n -1 )) (n -> 00); 

where the asymptotics follows from Stirling's formula. However, the quasi-optimal 
condition number k (m) behaves much better than just being of order 0(n 1//2 ). In 
fact, a comparison with (5.1) yields, as n — » 00 

5 97 

r (n) ~ n{n), k («) = 1 + — + ^^ +0 {n^), 
which is very close to optimal indeed. 

Example 7.5. We consider the example of Figure l.c, that is, the entire function 

f(z) = e^\ 

By the positivity of the Taylor coefficients, the maximum modulus function is 
also given byM(r) =/ _1 . A short calculation yields an explicit formula for the 
quasi-optimal radius, 

r = W(n), 

with the Lambert W-function as introduced in §2.2.2. To get our hand on the 
corresponding condition number bound, we realize that n\a n is the n-th Bell 
number whose asymptotics is well studied in the literature. Flajolet and Sedgewick 
(2009, Prop. VIII. 3) prove (using the concept of H-admissibility that we will study 
in §11) 

e e T «-l e e">-\ 

a n ~ = (n — > 00). 

r^lTir ^ + 1)^0 rl^lnnijo + X) 

Hence, asymptotically, we obtain the condition number bound 

K (n) = — ~ J2nn(r +1) ~ ^Innlogn (n -> 00), (7.5) 

where we have used the asymptotic expansion (de Bruijn 1981, Eq. (2.4.3)) 

W(t) =logf-loglogf + o(^^^) (f^oo). ( 7 .6) 
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Even though (7.5) looks like a possible, though moderate, 0(n ' (log n) ' ) growth 
of the condition number, things turn out to be much better than this. For instance, 
n = 100 yields the excellent quasi-optimal condition number k o (100) = 1.013. In 
§11 we will explain the surprising effect that K (n) is close to one for any order n, 
see Corollary 11.3. 



8. Entire Functions of Perfectly Regular Growth 

8.1. Order and Type of Entire Functions. Since r^(n) — > 00, an explicit asymptotic 
description of the optimization (7.3) requires a detailed study of the growth of the 
maximum modulus function M(r) as r — » 00. A fruitful characterization is by the 
order and type of /; for the following see Markushevich (1977, Thms. II.9. 2-9.5). 
The order p of an entire function f is given by 

log log M(r) .„ , 

p = limsup 6 , 5 — — =j0. (8.1 

r-^oo 1 logf 

Note that polynomials have order p = 0. If < p < 00 (which means that f is 
transcendental), the type t of / is given by 

log M(r) .„ , 

t = lim sup —5 — ±1 J> 0. 8.2 

r— xxj 1 fP 

We call / to be of minimal type if T = 0, of normal type if < T < 00, and of 
maximal type if t = 00. Order and type can also be read off from the coefficients a n 
of the Taylor series; if / is of order p, then 

n log n .„ , 

p — limsup, — —5 — (8.3) 
r n^oo^log(l/|fl„|) 

if / is of order p and type T, then 

T= —limsup n\a n \P /n . (8.4) 

ep n^oo r 11 

To arrive at an explicit asymptotic formula for r (w) (see Theorem 8.4) we 
need to consider a somewhat stricter class of entire functions (Valiron 1949, p. 45), 
though: an entire transcendental function of order < p < 00 is called to be of 
perfectly regular growth if the limit 

log M(r) 

t = hm (8.5) 

exists and is positive and finite; / is then of normal type T. The following funda- 
mental theorem is extremely helpful for the purpose of identifying such functions; 
for a proof see Valiron (1949, p. 108). 

Theorem 8.1 (Wiman 1916, Valiron 1923). Let f be an entire transcendental function. 
If f is the solution of a holonomic 17 differential equation of order q, then f is of perfectly 
regular growth with a rational order p ^ 1/q. 



I 7Holonomic differential equations are homogeneous linear with polynomial coefficients. 
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Example 8.2. The generalized hypergeometric functions 

p F q (h b p ;c t c q ;z) = f) f l)n ' ' ' f p \ n - (-b jr -c k ? W ) (8.6) 

n=0 v c l )n ■ ■ ■ \Cq)n «• 
are known to be (Luke 1969, §§3-3/5-1) 

• entire transcendental if and only if p ^ q; 

• satisfying a holonomic differential equation of order max(p, cj + 1). 

Thus, by Theorem 8.1, \i p ^ q, these functions are entire transcendental of perfectly 
regular growth with a rational order p ^ 1/ (<j + 1). It is an easy exercise in dealing 
with Stirling's formula 18 to calculate from (8.3) and (8.4) the order and type of 
these functions: 

1 

P = ZTTTn ' r=q+l-p (p<f). (8-7) 
H 1" 1 P 

Many transcendental functions can be identified as a generalized hypergeometric 
function (see Luke 1969, §6.2); if this relation is of the form 

f(z) = ctz* ■ p F q (bi, . . . , b p ; c lr . . . , c q ; f> z v ) (a, p £ 0, ft e N , v e N) 

then / is also of perfectly regular growth and we easily obtain, using (8.7), that 
the order and type of / are given by 

P = — v — , T={ q + 1 -p m v(^-r). 

1 ' V 

With the exception of the Airy functions, all the functions in the first section of 
Table 2 can directly be dealt with this way; it suffices to demonstrate just one such 
example in detail: 

cosz = ofi(;2'~5 z2 ) 

has p = 0, q — 1, v — 2, and f> = —1/4; therefore p — x — 1. 

Example 8.3. The Airy functions Ai(z) and Bi(z) satisfy a holonomic differential 
equation of second order, 

y"(z)-zy(z) =0. 

By the theory of linear analytic differential equations (Hartman 1982, p. 70), be- 
cause the leading coefficient of this equation is 1, the Airy functions are entire 
transcendental. Thus, Theorem 8.1 tells us that the Airy functions are of perfectly 
regular growth with a rational order p ^ 1/2. The precise values of the order and 
type can be read off from the asymptotic expansions (Abramowitz and Stegun 1965, 
Eq. (10.4.59-65)) of the Airy functions as z — »■ 00, which imply 

M(r) = - 7 ^ T7i J r3/2 (l + 0(r- 3 / 2 )) (r -> 00), 

with c = 1/2 for Ai(z) and c = 1 for Bi(z). Hence, by (8.1) and (8.2), we get 

3 2 
p=f, T=|. 



l8 Stirling's formula implies, for — c £ INp, that log | (c)„ | = n log )! — )! + 0(log n) as n — > 00. 
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Table 2. Various growth characteristics of some entire transcendental 
functions; all the functions with normal type are of completely regular 
growth and, a fortiori, of perfectly regular growth. The column for r (n) 
gives the asymptotics as n — ► 00. The angle 9 is understood to be restricted 
to — k ^ 9 ^ n. For l/r(z) the limit given is meant to be the interval 



(lim inf k ( 


w),limsupK («)). 


For the 


^-series ( — 


z;^)oo we assume 


that 




< q < 1. 
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8.2. The Asymptotics of the Quasi-Optimal Radius. A short calculation shows 
that any entire function f with the maximum modulus function 

log M(r) = rrP (p,r > 0) 

would have the quasi-optimal radius 

( n \ 1/p 

r«(n) = [-) . (8.8) 

By the definition (8.5), functions of perfectly regular growth satisfy the asymptotic 
relation 

log M(r) =rrP (l + o(l)) (r -> 00), (8.9) 

which suggests that (8.8) might still hold, at least asymptotically asn->oo. The 
following theorem shows that this is indeed the case; however, the proof is quite 
involved. 19 Concrete examples of the result can be found in Table 2. 



I9 Under the additional assumption of the non-negativity of the Taylor coefficients of /, it is possible 
to give a much shorter proof of this theorem; see Remark 12.2. 



ACCURACY AND STABILITY OF COMPUTING HIGH-ORDER DERIVATIVES 31 

Theorem 8.4. Let f be an entire transcendental function of perfectly regular growth 
having order p and type r. Then, the quasi-optimal radius satisfies 

( n\ 1/p 

r {n) ~ I — j (n — > 00). (8.10) 

Proof. The difficulty of the proof is to deal with the simultaneous limits r — >■ 00 
and n — >■ 00 whose coupling has yet to be established. To this end we introduce a 
transformed variable rj by 

We rewrite (8.9) in the form 

log(r""M(r)) = ne'?(l + o(l)) - -n - -log - = n •/»(»/)- - log -, 

p p X p T 

defining functions /n(j/) that satisfy 

=^(l+o(l)) -p-V 

note that the estimate o(l) holds locally uniform in // as n — » 00. By the properties 
of the maximum modulus function M stated in Theorem 7.1, we know that these 
functions /„ are strictly convex in rj and coercive, which means 

fn(v)^°° in -» ±°°)- 

The quasi-optimal radius r (w), which, by definition, minimizes r~"M(r), is now 
given in the form 

r«{n) = I — 1 , 

where rj n is the unique minimizer of /„ {rj). The assertion of the theorem is therefore 
equivalent to lim,,^^ n n = log p" 1 , which remains to be proven. 

Establishing the limit of n n proceeds by constructing a convex enclosure of f n 
for large n: for e > small, we define the strictly convex functions 

f ±£ (n)=e>1(l±e)-p- 1 n. 

The minimizer of f £ is explicitly given by 

1 

rj e = argmin/^/7) = log ( 1 + e ^ - 

Since f- e (t]) < f £ (rj) for all rj, and because /_ e is convex and coercive, there exist 
points fj and ^ e with r] £ < tj £ < rj £ satisfying 

f-e(%) =f-e(V e ) = fe(*le)- 

It is clear that rj ,rj e — > logp -1 as e -> 0; in particular, rj and 1J £ remain bounded. 
By the asymptotics of /„as)i4 00, we have, for n ^ n £ , 

fn(Ve) < feiVe) = f-e(%) ^ fn{%), 

fn(Ve)<fe(ye)=f-e(V e )^fn(V e ). 
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Thus, the strictly convex function f n is neither strictly increasing nor strictly 
decreasing between the points t] and fj £ . Hence, its minimizer t] n must lie there, 

V £ <Vn< Ve- 
Now, taking the limit n^co yields 

ri < liminf n„ -C limsupn,, < jT . 

Finally letting e — > proves that lim„^co = log p^ 1 as required. □ 

Remark 8.5. By means of (8.10) and (2.11) we can estimate the number of nodes m £ 
that a trapezoidal sum would need to achieve the relative approximation error e 
if we choose the quasi-optimal radius r = r (n). To this end we recall the Taylor 
series 

00 n 

W(z) = £ (-lf-V- 1 — (|z| < e- 1 ) (8.11) 

n=l n ' 

of the Lambert W-function (see de Bruijn 1981, §2.3) and obtain 

m £ f» en +jologe _1 . (8.12) 

Note how close this is already to the lower bound m > n given by the sampling 
condition (2.4). 

8.3. An Upper Bound of the Quasi-Optimal Condition Number. At a first sight 
the precise asymptotic description (8.10) of the quasi-optimal radius r («) does 
not tell us much about the size of the corresponding condition number K {n). In 
fact, restricting ourselves to subsequences of n which make the limes superior 
in (8.4) a proper limit, we just get 

logKo(n) ^ logKo(n) = o{n) (n — > 00). (8.13) 

Such a weak estimate could not even exclude a super-polynomial growth of the 
condition number. However, we can do much better (see the explicit asymptotic 
bound (8.18) below) by optimizing the upper bound 

M(r) 
\a n \r n 

from a dual point of view: by choosing the radius r in a way, such that the modulus 
of a n r n becomes maximal among all normalized Taylor coefficients; which directly 
leads us into studying the Wiman-Valiron theory of entire functions. For an 
account of the basics of this theory see Polya and Szego (1964, IV. 1-76); surveys of 
some more refined recent results can be found in Hayman (1974) and Gol'dberg, 
Levin and Ostrovskii (1997, Chap. 1.4). 

The fundamental quantities of the Wiman-Valiron theory are the maximum term 
of an entire function f with Taylor coefficients a n at a given radius r, defined by 

u(r) = max \a n \r", (8.14) 

n 

and the corresponding maximal index taking this value, called the central index, 

v{r) = max{n : \a n \r n = }i(r)}. (8.15) 

The asymptotic properties of these quantities are described in the following theo- 
rem; for a proof see Polya and Szego (1964, IV.68). 
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Theorem 8.6 (Wiman 1914). If the entire function f is of perfectly regular growth with 
order p and type r, then 

log M(r) ~ log^(r) ~ rr p , v(r) ~ Tpr p (r — > 00). 

We restrict ourselves to those entire functions / of perfectly regular growth for 
which eventually, if n is only large enough, each term \a„\r n (with a n 7^ 0) can be 
made the unique maximum term for a properly chosen radius. All the functions of 
Table 2 belong to this class. 

Remark 8.7. If a n 7^ for n large enough, then this property is known (see Polya 
and Szego 1964, IV.43) to be equivalent to the fact that \a n / a n+ ]\ becomes eventu- 
ally a strictly increasing sequence. This criterion is, for instance, satisfied by the 
generalized hypergeometric functions (8.6) with p ^ q: we find 



a n+l 



= (n + l) 



(n + ci) ■■■ (n + c q ) 



(n + b x )--- (n + b p ) 
which is therefore strictly increasing if n is only large enough. 

Thus, if a„ 7^ and n is large enough, then there will be a radius r n with 

n = v(? n ), \a n \?n = }i(r„). 

Theorem 8.6 yields the asymptotics (where n runs only through those indices with 

a n ^ 0) 

n = v(jn) ~ tpr n (n — > 00), 
which implies, in view of Theorem 8.4, the remarkable asymptotic duality 



n 



r n ~ ( ^- ) ~ u(n) (n -> 00). (8.16) 



We thus expect the bound (recall that r (n) is defined as the minimizer of K(n,r)) 

fc (n) _ %j < M(f H ) _ M(f n ) 

° v ; |««|r«(n)» " \a n \f n u(f n ) 

to be quite sharp for large n. Now, one of the deep results of the Wiman- Valiron 
theory is the following explicit bound of the ratio M(r)/f/(r) in general; for a 
proof see Hayman (1974, Thm. 6). 

Theorem 8.8 (Wiman 1914, Valiron 1920). Let f be an entire function of finite order p. 
Then, for each e > 0, there is an exceptional set E e of relative logarithmic density smaller 
than 1/ (1 + e) such that 



M{r) < p(l + e)}i(r)^27tlog]i(r) (r ft E e ). 

Shchuchinskaya (1982) has characterized those entire functions of finite order 
for which there are no exceptional radii, that is, for which E £ — 0. However, we 
did not bother to check her complicated conditions for any concrete functions. Let 
us simply assume the weaker condition that the sequence r n does eventually not 
belong to E £ for all e > 0. We would then obtain from Theorems 8.6 and 8.8, and 
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Table 3. For/(z) = Ai(z), a comparison of the quasi-optimal radius )"<>(«) 
with its asymptotic value (8.10) as taken from Table 2. This asymptotic 
value is already quite accurate for small n. The value of = \z n \ 

was actually computed by numerically solving the saddle point equation 
z«/'(zh)//(z«) = n in the complex plane. Note that lim„^oo Ko(n) = 
2/ v / 3 = 1.15470, see (10.15). 



n 


r*(n) 


K©(n) 


n 2/3 


K(«,n 2/3 ) 


1 


1.21575 


1.37413 


1.00000 


1.56499 


10 


4.72421 


1.19188 


4.64159 


1.21120 


100 


21.58047 


1.15832 


21.54435 


1.16003 


1000 


100.01668 


1.15506 


100.00000 


1.15523 



from (8.16) and (8.17), the asymptotic bound (where n runs only through those 
indices with a n 7^ 0) 

Ko(n) ^ K<>(n) ^ p^Jln log }i(r n ) ~ ^Inpn (n -> 00). (8.18) 

Note that this bound is consistent with the results obtained in Example 7.4 for 
/(z) = e 2 , in which particular case the bound of ic<>(tt) is even sharp; quite a 
success for such a general approach. In preparation of §10, we rephrase (8.18) by 
introducing yet another growth characteristics of /, namely the quantity 

< co = limsup K p^l < 1. (8.19) 
■ o:a„^o y/2npn 



n— 'rcau 



See Table 2, and also Polya and Szego (1964, IV.50), for some examples of 00. 

9. Relation to the Saddle-Point Method 

The results of the last section have shown that, for a certain class of entire 
functions of perfectly regular growth, the quasi-optimal condition number K (ft) 
grows at worst like 

1 ^ K<>(n) ^ K (n) = 0{n in ) (n -> 00). 

However, as we have seen in Examples 7.4 and 7.5, there are cases where the 
quasi-optimal condition number is asymptotically optimal, actually satisfying the 
best of all possible asymptotic bounds, K (n) ~ 1. We now develop a methodology 
which can be used to understand and prove this highly welcome effect for a large 
class of entire functions; concrete such examples will follow in the next sections. 

9.1. The Saddle-Point Equation. The key lies in the observation (Hayman 1974, 
Lemma 6) that the maximum modulus function M of an entire function / satisfies, 
except for a set of isolated radii (see also Theorem 7.1), the equation 



r— log M(r) 
ar 



r ,. = z Tz losfiz] 



z=z 

where zq G C is one of the points for which |zo| = tq and |/(zo)| = M(tq). We 
apply this observation to the quasi-optimal radius r n = r (n) which, by definition, 
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1/(7) z~ n \ 




1/(2) 2"" I 



Im( 2 ) 



a. f(z)=e>- 




Im(z) 



Figure 5. Plots of |z~"/(z)| for n = 31; left: /(z) = e 2 ; right: /(z) = 
Ai(z). The solid curve (red) is the image of the circle |z| = r («); showing 
that the maximum modulus along this circle is taken right at some saddle 
points. Note that the circle leaves these saddle points approximately in the 
direction of steepest descent. The left plot explains nicely the qualitative 
differences between the plots in Figure 4: where a circle gets close to 
being a level line of \z~ n f (z) , there must be oscillations of the integrand 
of the Cauchy integral. 



minimizes r n M(r). If not accidentally one of those isolated exceptions, this radius 
must fulfill the differential optimality condition 



r— log M(r) 
ar 



Thus, there is a complex number z n with 

r n = \z„\, M(r n ) = |/(z„)|, 



that satisfies the transcendental equation 



* S log/(z) 



l f(Zn)- 



This equation can be rewritten in the form 

F'(z n )=0, F(z)=z~ n f(z). 



(9-i) 



For functions F(z) that are analytic in a neighborhood of a point z n (with F(z„) ^ 0) 
it is well known (see de Bruijn 1981, §5.2) that F'(z„) = holds if and only if the 
modulus |F(z)| forms a saddle at z = z n . Since, by construction, |F(z)| has a local 
maximum at the saddle point z n in the angular direction, it must thus show a local 
minimum in the radial direction there; see Figure 5 for an illustration. On the other 
hand, by the convexity properties of the maximum modulus function M stated in 
Theorem 7.1, any saddle point z„ of |F(z)| satisfying |/(z„)| = M(|z„|) such that 
the saddle is oriented this way (local minimum in the radial direction and local 
maximum in the angular direction) will give us in turn the unique quasi-optimal 
radius r<>(«) = r n = \z n \. We have thus proven the following theorem. 
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Theorem 9.1. Let f be an entire transcendental function and let z n 6 C be a solution of 
the saddle-point equation F'(z w ) = with F(z) = z~ n f(z), that is, 

n=z n J —^-. (9.2a) 
f{z„) 

Ifz n = r n e ie » satisfies \f(z„)\ = M{\z n \), d ee \F{r n e' e »)\ < 0, and d„\F (r n e w ")\ > 0, 
then we get the following representation of the quasi-optimal radius: 

r {n) = |z„|. (9.2b) 

On the other hand, ifr (n) is a point of differentiability of M(r), then there is a solution 
z„ of the saddle-point equation that satisfies these three conditions. 

Theorem 9.1 allows us the actual computation of r (n); see Tables 3/4 and §10.4 
for some examples. 

9.2. The Saddle-Point Method. Taking the quasi-optimal radius r = r (n) we 
write the Cauchy integral (1.2) in the form 

1 f 2 * 

a » = 

2n Jo 

with F(z) = z~ n f(z). If |/(z)|, and thus |F(z)|, is small for those z on the circle 
that are not close to the saddle points z n of Theorem 9.1, the integral localizes to 
the vicinity of these saddle points and we can estimate 



F(r e 



ti„:z„=r e' 
saddle point 

It is actually possible to estimate each of the integrals 

1 r ... 1 
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/ F(r <> e ie )d9 = — [ e lo S F M e ) d6 
Je^e„ 2tz Je^e,, 



by the Laplace method (see de Bruijn 1981, §5.7). To this end we expand the 
function log/(re !e ) with respect to the angular variable 9; for # — S> we calculate 
that 

log f(re ie ) = log/(z*) + ja(z*)(0 - 6,) - \b(z*)(6 - 0*) 2 + O(0 - 0*) 3 (9.3a) 
with z* = re' 6 * and the coefficients 

a(z)=z^, b(z)=z«'(z). (9.3b) 

By specifying as the expansion point z* a saddle point z„ = r e' e " as in Theorem 9.1 
we thus have a{z n ) = n and therefore 

lo S F( r<> e ie ) =logF(z n ) - lb(z n )(6 - 6 n ) 2 +O(0 - 6„) 3 ; 

hence, by taking real parts, 

log |F(r e ie )| =log|F(z„)| - \Reb{z n )(6 - 6 n ) 2 + 0(9 - 6 n f. 

In particular, if |F(z)| takes when moving along the circle a strict local maximum 
at the saddle point z n , we infer that necessarily 

ReKz„)>0. (94) 
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Thus, the Laplace method is applicable and gives, by "trading tails", 



JL f e logF(r oe ' e 

27T Je^e„ 



In Je^e n 
1 

~ 2/r 



>gF(z„)-5b(zn)(e-e„) 2 



e io g F(z„)-ifc(z„)e 2 d0 = F ( z ") (9 5) 



y/2nb(z„ 

Summarizing our results so far, we get the following estimate of the Taylor coeffi 
cient a n : 

F(z) 



an " ^ Jj> 

saddle point 

Correspondingly we estimate the mean modulus by 
1 f 2n 



(9.6) 



2nJ 1^1^ 



|F(r e" 



U„:z„=r a e 



1 

2/r 



E 

n :z„=r o e ,e " 
saddle point 



saddle point 



gRelogF^) dQ 



I In 



a:z=r e'" 
saddle point 



(9-7) 



and, therefore, the quasi-optimal condition number by 



K<>(n) = K(w,r ) 



2/T 



2n 



F(r«e u 



y 

, ^ w y/Reb(z) 

0:z=r o e v v ' 

saddle point 

^ wr 



(9.8) 



W:z=r e' L 
saddle point 

As we will see in the following sections, for some interesting classes of entire 
functions our reasoning can eventually be sharpened by replacing the somewhat 
vague "~"-signs of approximation with rigorous asymptotic equality as n — » oo. 
Moreover, the estimate (9.8) is actually quite precise even for small n as is typical 
for such asymptotic estimates of integrals; see §10.4 for an example. 

9.3. Steepest Descent. In general, there is not much to further conclude about the 
approximate values of k (w) from the estimate (9.8). Thus, to get to a result like 
K o{n) ~ 1 we need some additional structure: a look at the examples of Figure 5 
tells us that there the circle of radius r passes through the saddle points of |F(z)| 
approximately in the direction of steepest descent. In the next sections we will 
explain why this is the case for some larger classes of entire functions. 

From general facts about the method of steepest descent in asymptotic analysis 20 
we learn (see de Bruijn 1981, p. 84) that the circular contour through the saddle 



20 For a detailed exposition see de Bruijn (1981, Chap. 5), Miller (2006, Chap. 4), and Flajolet and 
Sedgewick (2009, Chap. VIII). Gil, Segura and Temme (2007, §5.5) explain how steepest descent contours 
are used as an analytic tool for obtaining numerically stable integral representations of certain special 
functions; a topic that is certainly closely related to the theme of this paper. 
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point z„ is approximately of steepest descent if and only if b(z n ) is approximately 
real, that is, if and only if 

lmfc(z„)«0. (9.9) 

Note that this implies that the integrand in (9.5) has approximately constant phase. 
In fact, geometrically it is straightforward to see that the circle is the contour 
of steepest descent if and only if the off-diagonal elements of the Hessian of 
G(r,6) = log |F(re ,e )| vanish; at a saddle point z n = r e' e " as in Theorem 9.1 we 
actually obtain 



hessG(r , 9 n ) 



Refr(z„)r 2 
-hx\b(z n )r- 1 



-Imfr(z„)r 1N 
-Refc(zn) , 



(9.10) 



Now, assume additionally that the circle of radius r = r (w) passes through just 
one saddle-point z n (this amounts for the case Cl = 1 in §10.3). Then, we infer from 
the condition number estimate (9.8) and the steepest descent condition (9.9) that 



\F(z n ) 



K (n) 



V^eb(z n ) 



F(z„) 



\/b(z n ) 




f lmb{z n ) 
\Reb(z n ) 



1. 



This line of reasoning thus explains why the best of all possible results, k («) ~ 1, 
actually may come into place even though the radius r = r<>(n) itself was first 
introduced by optimizing just the upper bound ic(n,r) of the condition number. 



10. Entire Functions of Completely Regular Growth 

10.1. The Indicator Function. The reasoning of §9.3 relies on the remarkable fact 
(observed in Figure 5) that for certain functions the circle passing through the rele- 
vant saddle points is approximately tangential to the contour of steepest descent. 
This could be understood if F(z) = z~ n f(z) happens to grow predominantly in a 
radial direction. A first hint that this is exactly the right picture is the existence of 
the Phragmen-Lindelof indicator function 

h(6) =limsupr"^log|/(re' e )| (10.1) 

for entire functions of finite order p and normal type T. We recall some of its 
properties; see Markushevich (1977, §11.45) or Levin (1980, §1.15/ 16) for proofs: 

• h(9) is 2/T-periodic; 

• h(6) is continuous and has a derivative except possibly on a countable set; 

• if < p < 1/2, then < h(6) < t; if p > 1/2, then -t < h(6) < t; 

• t = rnaxg h(6). 

As it was convenient in §8 to consider the functions of perfectly regular growth, for 
which the limes superior in the definition (8.2) of the type T becomes the proper 
limit (8.5), we do the same with the limes superior in the definition of the indicator 
function here: 
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An entire function of finite order p and normal type T is called to be of completely 
regular growth (Levin 1980, Chap. Ill) if 

h(9)= lim to S 1/(^)1, ( 10 . 2 ) 

r^oo:r^E rP 

uniformly in 8. Here, the exceptional set E is required to have relative linear 
density zero; it will obviously be related to the zeros of /. In fact, if there are no 
zeros of / in an open sector containing the ray of direction 9, then (10.2) holds in a 
closed subsector without the need of an exceptional set. An important result of 
Levin (1980, p. 142) states that if (10.2) holds just pointwise for 8 in a set that is 
dense in [—71,71], then / is already of completely regular growth. This criterion 
can be used to check that all of the functions in the first section of Table 2 are of 
completely regular growth with the indicator functions given there: one just has to 
look at the known asymptotic expansions of /(z) asz^oo within certain sectors 
of the complex plane, as they are found, e.g., in Abramowitz and Stegun (1965). 
It is also known that the statement of Theorem 8.1 extends to completely regular 
functions, see Miiller (1997, p. 747). 

As developed mainly by Pfluger and Levin in the 1930s, there is a deep relation 
between the angular density of zeros of a function / of completely regular growth 
and the properties of its indicator function h(8). The following characterization of 
a density of zero will be of importance to us (Levin 1980, p. 155): 



# zeros z ^ r of f in an open sector containing the ray at 8q 
lim — - = 

r— >oo r" 

h(8) is jO-trigonometric in the vicinity of 8q; (10.3) 

where a function of 8 is called p-trigonometric if it is of the form a sin(p6 + ft) for 
some real a and f>. 



10.2. Circles Are Contours of Asymptotic Steepest Descent. We now look at a 
direction 0* in which there is the predominantly growth of /, that is, h(6*) = T. 
If there are at most finitely many zeros of / in an open sector containing the 
ray at 0* (which is the case for all of the functions in the first section of Table 2), 
then / will also be of perfectly regular growth and the indicator will be, by (10.3), 
p-trigonometric in the vicinity of 0* . In particular, we get 

/z(0*) = t, h'(d,) = 0, fe"(0.) = -t P 2 . (10.4) 

By the reasoning of §9 there will be a sequence z n = r e ie " (writing r = ?«(n) for 
brevity) satisfying the saddle-point equation (9.2a) with 8 n — >■ 0* as n — >■ 00. To 
show that the circle passing through z n is asymptotically a contour of steepest 
descent there, we look at the Hessian of log |F(z„)|. From (10.2) we first get 

log|F(z„)| =log|/(r e /e »)| -nlogr 

~ r^h{8 n ) — nlogr (n —> 00). (10.5) 
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Next, by Theorem 8.4 and (10.4), the Hessian of the right hand side, G(r,8) = 
rPh(9) — nlogr, becomes asymptotically diagonal: 

f r -2 \ 

hessG(r o ,0„) ~ np I J (n -> 00); 

note that this form of the Hessian is actually consistent with (9.10) and (9.4). Since 
the off-diagonal terms are zero, the 0-direction is, asymptotically, the direction of 
steepest descent. 

10.3. Condition Number Bounds. We follow the strategy of §9.2 and apply the 
Laplace method to the contour integral with radius r = r (n). However, instead 
of using the Taylor expansion (9.3) to simplify log F(re' e ) we now proceed by 
first recalling from §9.3 that contours of steepest descent yield integrands of an 
asymptotically constant phase and by next using the indicator function (10.2) 
to simplify log |F(re ,9 )|, asymptotically as r — >■ 00. Note that the Laplace method 
rigorously applies if there is a proper decay of log \ F(re lB ) |, as r — » 00, for directions 
8 far off those that belong to the saddle points. Assuming this to be the case for 
the given / (it can be checked to be true for all the functions in the first section of 
Table 2), we get for the Cauchy integral (1.2), because of (10.5), (10.4) and (8.10), as 
n — > 00: 



2nr. 



TYq J — tt 27T J — 71 

!_ y e ilmlo S F(r e w ) e RelogF(r oe w ) + y o (t-e) 2 h''(8) dt 

-7Z a , J —00 
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- E e- me f{r«e ie ). 



r <> d:h{6)=T 



Likewise, we get, as n — > 00, 



and certainly 



27T?-" j-tt 
1 

^J2npn ■ 

M(r ) 



[ n \f{rj e )\de= ^- ( n e Rel °s F ^ eie Ue 

J — 71 2TL J —71 



r o 6:h(d)=T 



~r-« max \f{rj e )\. 
n e-.h{e)=r 

To summarize, we have proven the following theorem. 

Theorem 10.1. Let f be an entire function of completely regular growth with order p, 
type t, and Phragmen-Lindelof indicator function h(8). If f has at most finitely many 
zeros in some sectorial neighborhoods of those rays of direction 6 for which h{8) = r and 
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if I/) decays properly, for large radius r, in the angular direction off these rays, then we 
have 

max 1/(^)1 



and 



Ko(tt) 

yjlnpn 



Ko{n) 



E e 

8:h(9)=T 

E 

e:h(e)=r 



ine f(r»e ie ) 
\f(rJ d )\ 



(n 



a n £ 0) 



(10.6) 



(n — > 00 : a n ^ 0). 



(10.7) 



£ e- me f(ue ie ) 

8:h(6)=T 

That is, the quasi-optimal condition number K (n) of the Cauchy integral is asymptotically 
equal to the condition number of the finite sum YLe-.h(6)=r e~ in6 fir^e 1 ®). 

Let us introduce the number of global maxima of the indicator function, 

Q = #{0 : -7T < 6 < n,h(e) = r}. (10.8) 

Now, by Theorem 10.1, Q = 1 clearly implies that lim„^co K (n) = 1 and that the 
quantity defined in (8.19) satisfies co = 1; this observation is precisely matched 
by two examples in Table 2. On the other hand, if fi > 1 then it seems, at a first 
sight, that the condition number of the finite sum Y,e-h(e)=T e ~' ne f(rc,e' e ) could 
suffer from severe cancelation. However, as the next theorem shows, there will be 
generally no such cancelation for the class of functions considered in this section. 
(But see §10.4 for an example of severe resonant cancelations in a different setting.) 

Theorem 10.2. Let f be an entire function of completely regular growth which satisfies 
the assumptions of Theorem 10.1 as well as those that led to (8.19), that is, to co ^ 1. Then, 
this bound can be supplemented by 



< fT 1 ^ liminf 



^ limsup 

:«n/0 y/2npn n-i>oo:fl„/o y/2npn 



co < 1, 



(10.9) 



and the quasi-optimal condition number K (n) is asymptotically bounded as follows: 

1 ^ liminf K (n) ^ limsup K (n) ^ Q • co. (10.10) 



In particular, we have 

cr 1 



CO 



v 1 \ v K o( n ) 

lim Ko(n) = lim ——= 

n->oo:« n5 ^o n->oo:fl„^o Ulnpn 



1. 



Proof. The obvious estimate 

~ ine f{rJ e ) 



E * 

6:h(6)=T 



(lO.ll) 



< E l/(r/)l<n- max \f{r^)\ 
e-.h(e)=r e-.h(e)=T 

yields, by Theorem 10.1 and (8.19), the asymptotic bounds asserted in (10.9). 
Moreover, (8.19) and (10.6) imply 



lim sup 



max \f(r <> e l[ 
0:h(6)=r 



E * 

8:h{6)=T 



-in0 



f(r«e U 



= 00 1. 
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Hence, by using (10.11) once more to estimate the numerator in (10.7), we get 



J0\ 



max f(r e 

, . , 8:h{0)=T 

limsup Ko(n) 12 limsup — — — ^ 12 ■ co, 



n— >oo:a„^0 n— >oo:a n ^0 



£ e- ine f(rj e ) 

6:h(e)=T 



which proves the asserted asymptotic bound (10.10). □ 

Example 10.3. If, by the symmetries of the function f in the complex plane, there 
is just one single phase (p„ G 1R that allows us the representation 

e i( f" ■e- ine f{rj e )= max \f {rj e *)\ (10.12) 

8 t :h(0*)=T 

for all 6 with /z(0) = t, then we get by Theorem 10.1 that already the best of all 
possible bounds holds, namely 

lim K (n) = 1, lim * o( w ) _ q-1 (10.13) 

We than have, by definition, a; = Q . Note that the symmetry relation (10.12) 
applies to all of the functions of the first section of Table 2, except for the Airy 
functions Ai(z) and Bi(z) which will be dealt with in the next two examples. 

Example 10.4. The point of departure for discussing the Airy function Ai(z) is 
the asymptotic expansion (Abramowitz and Stegun 1965, Eq. (10.4.59)) 

. . , \ 1 _ 1/4 2 2 3/2 S (-l) k T(3k + \) _ 3k/2 . 

Al(z) ~ — z 1/4 e 3 Z > - — r— -^-z (z 00 : arez < n). 

1 ' 2n £5 9 k T{2k + l) v 

(10.14) 

This implies, by Levin's criterion given above, that Ai is of completely regular 
growth. Moreover, we get 

|Ai(re ,e )| = ^r- 1/ V3 r3/2cos (2 e )(l+0(r- 3/2 )) (r ->■ co : |0| < tt), 

from which we can directly read off the order p — 3 / 2, the type T = 2/3, and the 
Phragmen-Lindelof indicator function 

/i(0) =-§cos(§0) (|0|<tt). 

Note that this indicator 7z(0), continued as a 27i-periodic function, is p-trigonometric 
exactly for 6 7^ krc (k G Z). Thus, by Levin's general theory, there is a positive 
density of zeros in an arbitrary small sectorial neighborhood of the ray at 8 = — n; 
indeed, Ai(z) has countably many zeros along the negative real axis and no zeros 
elsewhere. We have h(6) = r for 6 = ~^-\ n ) hence (1 = 2. The expansion (10.14) 
implies for these maximizing angles that 

Ai(re ± i 7ri ) = ^Lr- 1/4 e3 r3/2 (l + 0(r~ 3/2 )) (r co), 
2v 7t 

that is 

argAi(nr t § ,ri ) =^7 + 0(r~ 3/2 ) (r -> co). 
6 
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Hence we obtain, because of h(— j7t) = h(^n): as n — > oo, 



-in8 



Ai( r<> e w ) 



8:h(6)=T 



2 . 71 . 2 . 71 

„^TTni„-rl | — 57771! — -ri 

e3 e 6 + e 3 e 6 



Ai r e3 7 



2|cos(f + |tto) 



Ai r e3 7 



and 



E \ Ai (r«e U 



e-.h(8)=T 

Now, 



Ai r e3 7 



max |Ai(r e /e )| 

9:h(8)=T 



Ai r e3 ; 



|cOs(f + §7TO) 



'V3/2 




n ^ 2 (mod 3), 
n = 2 (mod 3), 



in accordance with the fact that the Taylor coefficients of Ai(z) satisfy a n 7^ if 
and only if n ^ 2 (mod 3). Altogether, Theorem 10.1 gives us then 



lim jc (n) = — =, 

«^oo:fl„^0 V3 



a> = lim —== = —=. 

fWoo:a„^0 y'lnpn V3 



(10.15) 



We observe that the general upper bound given in (10.10) is sharp here. An 
illustration of the limit result (10.15) by some actual numerical data for various n 
can be found in Table 3. 

Example 10.5. As for Ai(z) in the last example, the discussion of Bi(z) begins 
with its asymptotic expansions (Abramowitz and Stegun 1965, Eq. (10.4.63-65)) 
as z — » 00 in different sectors of the complex plane. Skipping the details, we 
get that Bi is of completely regular growth with order p — |, type T = |, and 
Phragmen-Lindelof indicator 



h(e) = l\cosm\ 



< n). 



Thus, h(6) = t for 6 = ±§7r and also for 9 = 0; hence Q = 3. The asymptotic 
expansions yield 



Bifre^ 7 "') 



e 3 
— =7 
2^ 



A r 3/2 



e3 r ~"{l+0(r- 3/2 )) (r^co) 



and 



Bi(r) 



-l/4j^ 2 (l + 0(r -3 2 



7T 



(r" 3/z )) (r ->«>), 



that is arg B(r) = and 



argBi(re ± 3 m ') = ±f +0(r~ a/z ) 



3/2^ 



00 . 



Hence, as r — >• 00, 



|Bi(re d 



i|Bi(r)| 
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Table 4. For /(z) = Bi(z), a comparison of the quasi-optimal radius )"<>(«) 
with its asymptotic value (8.10) as taken from Table 2. This asymptotic 
value is already quite accurate for small n. The value of = \z n \ 

was actually computed by numerically solving the saddle point equation 
Znf (z n ) / f(z„) — n in the complex plane. Note that Iim n ->oo K<>(n) 
4/3 = 1.33333, see (10.16). 



n r<>(w) Ko(n) tt 2 K(n,n 2 ' 3 ) 



1 


1.36603 


1.35408 


1.00000 


1.57640 


10 


4.72421 


1.37605 


4.64159 


1.39833 


100 


21.58047 


1.33751 


21.54435 


1.33948 


1000 


100.01668 


1.33375 


100.00000 


1.33394 



and thus, as n — >■ 00, 



£ e-" ,e Bi(r e'' e ) 

9:/i(0)=t 



1 2 7T. 1 2 . 7T. 

— e3 e 3 + -e 3 e 3 + 1 



|Bi(r 



and 



Now, 



E l Bi (^' 

8:h(e)=T 



= |l + cos(f - §7rn)| ■ |Bi(r )|, 
2|Bi(r )|, max |Bi(r e l9 )| ~ |Bi(r 

6:h{0)=T 



|l +COs(^ — J 7TW) I 



3/2 n ^ 2 (mod 3), 

^0 n = 2 (mod 3), 

in accordance with the fact that the Taylor coefficients of Bi(z) satisfy a n 7^ if 
and only if n ^ 2 (mod 3). Altogether, Theorem 10.1 gives us then 

4 jc(n) 2 



lim Ko(n) 



a; = lim = -. (10.16) 

n->oo:a„J=0 yjlnpn 3 



An illustration of the limit result (10.16) by some actual numerical data for various 
n can be found in Table 4. 



10.4. A Resonant Case: /(z) = 1/T(z). In the statement of Theorem 10.1 the 
condition on the zeros of / cannot be disposed of: if / possesses infinitely many 
zeros in the vicinity of its directions of predominant growth, then it may happen 
that a pair of saddle points recombines in the limit r — » 00 to a single maximum of 
the indicator function h{&). That is, even though we have CI = 1 in the limit, the 
contributions of the two saddle points may yield resonances in (9.8) as n — » 00; 
thus K (n), as well as K*(n), may behave quite irregular. 

We demonstrate such a behavior for the entire function f(z) — 1 /T(z), whose 
zeros are located at 0, —1, —2, —3, . . . This function has order p — 1, but is of 
maximal type T = 00 (see Levin 1980, p. 27). Therefore, at a first sight, the results 
so far do not seem to be applicable at all. However, using Valiron's concept of a 
proximate order p(r) it is possible to extend the definition of functions of completely 
regular growth and of their indicator functions in such a way that the results cited 
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Kjn) density F'(t) 




500 1000 1500 2000 2500 -12 -10 -8 -6 -4 -2 



a. quasi-optimal condition number K<,{n) (1 ^ n ^ 2600) b. histogram of log log Ko(«) 

Figure 6. Left: plot of the quasi-optimal condition number K (n) (1 ^ 
n ^ 2600); /(z) = l/r(z). Within the shown range of n, the maxi- 
mum is taken for n = 2006: k o (2006) = 47067.2. Note that there is 
not much to be gained from using the optimal radius r*(n) instead 
of u{n): k*(2006) = 47063.9. Right: plot of the density histograms of 
t = loglogKo(M) (1 < n < N) for N = 100000 and N = 1000000 and of 
the density F'(t) belonging to the distribution F(f) = ^ arccos(exp(— e 1 )), 
printed transparently on top of each other. Since there is such a close 
agreement we are led to conjecture the limit law (10.20) and, therefore, 
lim inf„_i.oo K«(n) = 1 and limsup H _i.oo Ko(ft) = °°- 



above still hold true (see Levin 1980, §1.12). By Stirling's formula, and Euler's 
reflection formula 

r(z).r(i-z) = ^- T , 

sm(7rz) 

we get the following asymptotic expansion, valid uniformly in 9: 

log|l/r(re /e )| „ cos0 + 0sin0 nt . . r . . 

& ' v '- 1 = -cos6-\ : hO r 1 ) (r ->■ 00 : r 4. £), (10.17) 

r log r log r 

where the set E of possible exceptions has relative linear density zero. From this 
we can read off that 1/T(z) is a function of completely regular growth with a 
proximate order p(r) given by rP^ = r logr; the indicator function is then 

h(6) = -cos0. 

Now, the problem is that this indicator becomes asymptotically maximal at the 
single direction 8 = ±7r, which is actually the direction of the ray that contains 
the countable many zeros of 1/T(z). In fact, a closer look at (10.17) reveals that 
this single maximum is formed, in the limit r->oo, through a recombination of 
two distinct maxima for finite r. And indeed, Figure 6.a shows quite an irregular 
behavior of the quasi-optimal condition number K (n) (the picture would be 
essentially the same for the optimal condition number K*(n) itself, though much 
more difficult to compute). 

The quasi-optimal radius r (n) can straightforwardly be obtained by means of 
the saddle-point equation (9.2a): that is, J-<>(n) = \z n \ where z n is one of the two 
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complex conjugate solutions of 

n = 4 1 ° § f(l) = - 2!/ ' (z); 
we choose Imz„ > for definiteness. Asymptotically, as n — »■ oo, this saddle- 
point equation can actually be solved explicitly in terms of the principal branch 
of the Lambert W-function: using the asymptotic expansion (Abramowitz and 
Stegun 1965, Eq. (6.3.18)) of the digamma function ip we obtain 

—zip{z) = -zlogz+ - + 0(z~ 1 ) (|argz| < n), 
and therefore, as n — »■ 00, 

W(i-n) 

which we take as the definition of the radius r n and the angle 71/ 2 < 8 n < n. 

A detailed quantitative analysis of K (n) can now be based on the well-known 
fact (see Hayman 1956, p. 91) that the saddle-point analysis of §9.2 is applicable to 
f(z) = 1/T(z): in fact the approximations (9.6) and (9.7) are asymptotic equalities 
as n — > 00. We find that they can be recast in the form 

T~ \1/T(r n e w ") 
nn r" 



■ COS <p n , (10.19a) 



_ , . / 7ZYI 1 . . , , 

Ko(n) ~ W — -|sec0„|, (10.19b) 

K<>(n) ~ I sec cp n I (10.19c) 
with the collective phase approximation 21 

<pn - (n - ^ 6 " -0 n + 12 ^" 1^2 ^ ~ ^arccot (cot0 M - 6„ esc 2 6 n 

The asymptotics (10.19c) does not only explain the very possibility of resonances, 
it actually gives excellent numerical predictions even for rather small values of n 
such as those illustrated in Table 5. 

Based on Table 5 and Figure 6. a it is certainly quite reasonable to conjecture that 
liminf„_ 5l00 K (n) = 1. On the other hand, by just looking at the rather randomly 
distributed positions n of the resonances and the corresponding extreme values of 
K^(n) we could not really establish any serious conjecture about the probable value 
of limsup„_i>oo K (n). Instead, we look at the statistics of the values of k (w) for 
1 ^ n ^ N. The very close agreement of the two histograms shown in Figure 6.b 
suggests that there should be a limit law of the form 

lim AT 1 -#{1 ^ n ^ N : log log K (n) < t} = F(t). (10.20a) 

N— >co 



2I Hayman (1956, p. 92) basically states the same results with the much simpler phase approximation 

<p„ = {n-Die^sW-en-Bn), 

which is, however, numerically far less accurate for small values of n and would not allow such a 
precise prediction of K (n) as in Table 5. 
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Table 5. The precision of the asymptotics (10.19c) near some resonances. 



11 


i^O\ n ) 


| SCL^fl J \ 


11 




Cpf I d)*.\\ 
| SCL ^ If/ft j | 


2002 


1.018 


1.018 


10931 


1.006 


1.006 


2003 


1.034 


1.033 


10932 


1.124 


1.124 


2004 


1.301 


1.300 


10933 


1.498 


1.497 


2005 


2.354 


2.352 


10934 


2.798 


2.797 


2006 


47067.162 


42811.637 


10935 


138149.749 


143720.416 


2007 


2.355 


2.353 


10936 


2.798 


2.797 


2008 


1.301 


1.300 


10937 


1.498 


1.497 


2009 


1.034 


1.033 


10938 


1.124 


1.124 


2010 


1.018 


1.017 


10939 


1.006 


1.006 



If the phases <p n were equidistributed modulo n (and the empirical data of the first 
one million instances strongly point into that direction) we would immediately 
find from (10.19c) that the distribution would be 

F(f) = — arccos(e~ e ). (10.20b) 

In fact, we observe that the thus given density F'(t) is very well approximated by 
the histograms in Figure 6.b and we therefore conjecture that the limit law (10.20) 
is correct. Now, since F'(t) > for all f 6 K, this conjecture would also imply that 

liminf K (n) = 1, limsupKo(n) = 00. 

n— >oo n— ^oo 1 

Actually, things are not as bad as such a spread of the condition number might 
suggest: from ~ arccos(l/100) = 0.9936 we infer that just about 0.64% of all n (in 
the sense of natural density) have k («) 100; that is, as much as at least 99.36% 
of all the Taylor coefficients a n enjoy to be computed with a loss of less than two 
digits. We find that the asymptotic median of k (w) would be as small as y/2. 

Remark 10.6. In the same vein, a worst-case analysis based on Figure 6. a tells us 
that there will be just a loss of at most three digits in computing the first one 
thousand of the Taylor coefficients of 

by means of a Cauchy integral with radius r<>(n). Note that the only competitor of 
this approach, namely using the recursion formula (see Luke 1969, §2.10) 

n-1 

a = 0, a\-l, a 2 = 7, a n = na 1 a n - fl 2 «n-l + {- 1 ) k ^) a n-k (» > 2), 

k=2 

is much worse behaved and suffers from severe numerical instability almost right 
from the beginning: in hardware arithmetic all the digits are lost for n ^ 27. 
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11. H- Admissible Entire Functions 

The function f(z) = e eZ_1 of Example 7.5 is not covered by our results so far: it 
has order p — 00. Nevertheless, the general idea of using the saddle-point method 
(see §9) can certainly also be applied to functions that grow even stronger than /. 
Hayman (1956) has axiomatized an important class of functions (with predomi- 
nant growth in the direction of the real axis), for which the saddle-point method 
is applicable along circular contours and which enjoys nice closure properties. 
Expositions of this method can be found in Wong (1989, §11.79), Odlyzko (1995, 
§12.2), Wilf (2006, §54), and Flajolet and Sedgewick (2009, §VIII.5). 

Hayman's method is based on the Taylor expansion (9.3) with the expansion 
point z* = r, that is, on the Taylor expansion 

log f{re ie ) = log/(r) + ia(r)6- ^b{r)6 2 + O{0 3 ) (0 0), (11.1a) 

where the coefficients are given by 

f(r) d 

air) = r . = r— log fir), bir)—ra'ir). (11.1b) 
f{r) dr 

Now, an entire function f(z) that is positive on (tq, 00) for some r$ > is said to 
be H-admissible, if it satisfies the following three conditions: 

• b(r) — >■ 00 as r — >■ 00; 

• for some function #o( r ) defined over (tq, 00) and satisfying < #o( r ) < 71 1 
one has, uniformly in |0| ^ $o( r )/ 

• uniformly in 8q (r) ^ \8\ ^ n 

fire-) = °Jfg±. 

However, one rarely checks these conditions directly but relies on the following 
closure properties instead. 

Theorem 11.1 (Hayman 1956). Let f and g be H-admissible entire functions and let p 
be a polynomial with real coefficients. Then: 

(a) the product f(z)g(z) and the exponential ef^ are admissible; 

(b) the sum f(z)+p(z) is admissible; 

(c) if the leading coefficient of p is positive then f(z)p(z) and p(f(z)) are admissible; 

(d) if the Taylor coefficients ofePW are eventually positive then e v ^ is admissible. 

For instance, with the help of this theorem it is fairly obvious to see that the 
functions fiz) = e z and fiz) = e eZ_1 are both H-admissible. On the other hand, 
the H-admissibility of functions like fiz) = z _i ^ 2 J; c (2 v / z) has to be inferred more 
labor-intensive from the definition. 

From the definition of H-admissibility we immediately read off that the maxi- 
mum modulus function is given, for r large enough, by 

M(r)=f(r), (11.2) 
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which, by the strict convexity of log M(r) with respect to logr, 22 by Theorems 7.3 
and 9.1, implies that the quasi-optimal radius r = r (n) is the unique solution of 

a{r<>) = n (11.3) 

for n large enough. Hayman's main results are summarized in the following 
theorem. 

Theorem 11.2 (Hayman 1956). Let f be an entire H-admissible function. Then, for the 
quasi-optimal radius r = r (n), we have 23 

'-'-"^-^Mm ( ™ ); (IM) 

in particular, we get a n > for n large enough. Moreover, we have, uniformly in the 
integers n, 24 

a nr n 1 / / {n-a(r)f\ M \ . . , 

ex P ~ m\ (r^oo). (11.5) 



f(r) J2nb{r) V V Mr) 
Finally, the ratio a„/a n+ \ forms an eventually increasing sequence since 25 

r (n) ~ ~ [n — > 00). (n-8) 

As for the conditions numbers, we straightforwardly get the following corollary; 
for reasons of a better comparison we have also included the quantities of the 
Wiman-Valiron theory as introduced in §8.3 (their asymptotics can directly be 
read off from (11.5)). 

Corollary 11.3. Let f be an entire H-admissible function. Then 

v i \ 1 1 ■ 1 
lim Ko(n) = 1, lim — — = 1. 

n-^co n^>oo y /2nb(r <> (n)) 

Moreover we have v(r) ~ a{r) as r — > 00 and 

M(r<>(n)) ~ yj2nb(re(n)) }i(ro(n)) (n ->■ 00). 



22 Note that this strict convexity implies 6(r) = ('"^:) 2 log/(r) > for all r > r$. 

2 3Note that (11.4) can be thought of as being a generalization of Stirling's formula, cf. Examples 5.1 
and 7.4: this was the original headline of Hayman's (1956) work. 

2 4p>ecause of f(r) = EtLo a k rk t the quantities a„r" / f(r) form, if a H ^ for all n, a probability 
distribution in the discrete variable n. The result (11.5) thus tells us that this probability distribution is 
asymptotically, in the limit of large radius r — > oo, Gaussian with mean a(r) and variance b(r). 

2 5Note that the asymptotic representation (11.8) of the quasi-optimal radius holds for the generalized 
hyperbolic functions (8.6) with p ^ q, too: namely, we have by Theorem 8.4, Example 8.2 and Remark 8.7 
that 

r«(«)~n'J+ 1 -P~ — ^ (n -*•«,). (11.6) 

On the other hand, such a representation is not valid for the function of Example 12.4. However, there 
the following corollary of (11.6) is nevertheless correct: 



a n-\ 



(n — > 00). (11.7) 



Hence, if we restrict ourselves to those n for which 8 n _i, a n +\ / 0/ we observe that (11.7) does in fact 
hold for all the functions of Table 2 but the function 1/I*(z). Whether this fact is just a contingency or 
whether it is for some deeper structural reason, we do not yet know. 
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Applications have already been discussed in Examples 7.4 and 7.5. 

Remark 11.4. If the entire H-admissible function f is of finite order p with normal 
type t, it is instructive to compare Corollary 11.3 with Theorem 10.2. From the 
definition of H-admissibility it then follows that: 

• / is of perfectly and of completely regular growth; 

• there is just one direction of predominant growth, O = 1 with h(0) = t; 

• / has at most finitely many zeros in the vicinity the positive real axis. 

Thus, / satisfies the assumptions of Theorem 10.2 and also those that have led to 
the definition (8.19) of to. Therefore, by CI = 1 we get from (10.9) and (10.10) that 
co = 1 and 

lim K<>(n) = 1, lim = 1. 

(Recall that H-admissible functions have a„ > for n large enough.) Further, by 
Theorem 8.6 we have v(r) ~ rprP. These results are consistent with Corollary 11.3; 
a comparison gives, by using (8.10), the asymptotic equations 

a(r) ~ rpr p , b(r) ~ Tp 2 r p (r — > 00). (n-9) 

Formally, as suggested by (11.1b), these equations could have been obtained from 
differentiating the asymptotic equation log/(r) = logM(r) ~ rrP (which just 
states the perfectly regular growth of the function f, see (8.5) for the definition). 
The differentiability of these asymptotic equations has also been observed by Polya 
and Szego (1964, IV. 70/71) under the weaker assumption that f is a function of 
perfectly regular growth with a n ^0. 

12. Entire Functions with Non-Negative Taylor Coefficients 

In this final section we consider entire transcendental functions / which have 
non-negative Taylor coefficients: a n ^ for all n. Such functions are typically met 
as generating functions in combinatorial enumeration or in probability theory. The 
non-negativity of the Taylor coefficients implies at once that 

M(r)=f(r) (r>0). (12.1) 

Thus, by Theorem 7.1, we infer that log/(r) and hence log(r~"/(r)) are strictly 
convex functions of log r. Moreover, since 

f^r- n f(r) = r~ n - 2 £ (n + 1 - k) (n - k)a k r k > (r > 0), 
dr fc=o 

we conclude that the function r~ n f(r) itself is strictly convex, too. The same 
reasoning that led to (11.3) in the last section proves the following simplification 
of Theorem 9.1. 

Theorem 12.1. Let f be an entire transcendental function with non-negative Taylor 
coefficients: a n ^ Ofor all n. Then, the quasi-optimal radius r = r (n) is given as the 
unique solution of the convex optimization problem 

r = argmin r~"/(r), (12.2) 

r>0 
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and, equivalently, as the unique solution of the real saddle-point equation 

r ° = " (r ° > 0) ' (12 ' 3) 

Remark 12.2. If / is a function of perfectly regular growth (of order p and type t) 
with non-negative Taylor coefficients, Theorem 12.1 yields the assertion of Theo- 
rem 8.4 with a proof that is much shorter than the one of the general result given 
there: first, from the definition of perfectly regular growth in (8.5) and from (12.1) 
we get 

log/(r) ~ rr p (r — > 00); 

next, since the Taylor coefficients are non-negative, we may differentiate this 
asymptotic equation (see Polya and Szego 1964, IV. 70) and obtain 

f'(r) 

r r . . ~ Tpr p (r — > 00); 

fir) 

therefore, by recalling r (n) — » 00 as n — » 00 (see Theorem 4.6), the saddle-point 
equation (12.3) is asymptotically solved by 

n N l/ P 



r oi n ) ~ I — J ( n -> °°)r 

which is, finally, the assertion of Theorem 8.4. 

Example 12.3. We pick up the computation of the gap probabilities E2(n;s) as 
discussed in Example 3.1, this time striving for small relative errors instead of 
absolute errors. As we have seen, the generating function is given by the Fredholm 
determinant 

00 

/CO = E £ 2(A:;s)z i: = det(7-(l-z)X| L 2 (0/S) ), K(x,y) =wDc{n{x-y)), 
k=0 

which is known to be, as a function of z, an entire function of order p — O. 26 Fig- 
ure 7.a shows that, for n = 10, the quasi-optimal radius r (as computed from (12.2) 
by means of Matlab's f minbnd command) varies over about 20 orders of magni- 
tude as the parameter s runs through the interval 2 $C s $C 18. The corresponding 
condition number satisfies k = 1, up to machine precision throughout. We will 
explain this optimal condition number result and the strong variability of the 
radius by discussing a "mock-up" model in the next example. Note that even 
though Figure 7-b shows a significant accuracy improvement in the tails, we do 
not get the full accuracy that we would have expected from k = 1. The reason 
is simply that the numerical evaluation of the Fredholm determinants does not 
satisfy the model assumption (3.2); see Bornemann (2010, §4). 



2 Generally, if the kernel K(x,y) satisfies a Holder condition with exponent a, with respect to either 
x or y, then f(z) = det(7 — zK) is an entire function of order p ^ 2/(1 + 2a); see, e.g., Lax (2002, 
Lemma 24.10). 
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a. quasi-optimal radius r = r c as a function of s 



Figure 7. Left: the quasi-optimal radius r as a function of s for calcu- 
lating the gap probability E2(10;s) of GUE as the 10-th Taylor coefficient 
of a Fredholm determinant. Right: the relative error of the calculation; 
the dotted line (red) shows the errors for the radius r = 1, the solid 
line (blue) shows the errors for the radius r (see also Example 3.1 and 
Figure 2). Note that, though k = 1 throughout the range of s, there is 
still a noticeable loss of accuracy in the tails: this is because the Fredholm 
determinant is not computed to small relative but small absolute errors; 
hence the model assumption (3.2) is violated. Nevertheless, r<> gives a 
significant improvement over the fixed radius r = 1 that belongs to the 
concept of absolute errors. 



8 10 12 14 16 18 

S 

b. relative error 



Example 12.4. Lacking a proof that k ~ 1 in Example 12.3 we analyze a "mock-up' 
Fredholm determinant in full detail, namely the ^-series 



naW) (o<?<i). 

k=0 



/(z) = {-z;q) c 

By a result of Euler (see Andrews et al. 1999, Cor. 10.2.2) it is known that 

m 



f 1® 



where 

(q; q ) k =(l- q )(l- q 2 )---(l-q k ); 

in particular, / has positive Taylor coefficients. By using (8.3), / is easily seen to 
be an entire function of order p = 0. Now, a numerical experiment shows that 
7Co(n) = 1 up to machine precision for n = 20 and q = 1/2. Hence we aim at 
proving that K (n) — >■ 1 as n — > 00. 

A natural first try would be to check / for H-admissibility, see Corollary 11.3. 
This approach is doomed to fail, however, since we get the following asymptotics 
from an application of the Euler-Maclaurin sum formula: 



a(r) 



fir) 
f(r) 



St 



k=0 



logr 1 

V " log(l/qr) + 2 



0(r- 



(12.4) 



ACCURACY AND STABILITY OF COMPUTING HIGH-ORDER DERIVATIVES 



53 



but 



which remains bounded (recall that H-admissibility would require b(r) — > oo). 
Nevertheless, from (12.3) and (12.4) we get the strong estimate 

r»(n) =q 1/2 -" + 0{l) (n -►«>), (12.5) 

which not only shows that r (n) grows exponentially fast for this slowly growing 
function f(z), but also that r varies very strongly with respect to the parameter q 
(an effect that we had already observed in Example 12.3). 

To study Ko(n) we look more closely at log/(re' e ) for large values of r. Applying 
the Euler-Maclaurin sum formula once more and using a uniformity criterion of 
Levin (1980, p. 142), we get that 

„ , ,-/ 1 loe(r) 2 1, 1 , .„ , . 7r 2 — 39 2 „. 1s 

Relog/(re 19 ) = - *\ > + -logr + - log(l/f) + , Wl . , +0(0, 



21og(l/<?) 2 fa '12 m " 61og(l/«/) 

and 



(12.6a) 



Imlog/(r e ,e ) = Q - Xo f" - + + ^- S in(0)r- 1 +0(r- 2 ), (12.6b) 
log(l/^) 2 1 — q 

uniformly in 6 as r — >■ 00 (r ^ E) with the possible exception of a set E of relative 
linear density zero. The first asymptotics, (12.6a), means that / is of completely 
regular growth with the proximate order p(r) (see Levin 1980, §I.i2), 27 

Mr) = 1 Jog(r)i 

21og(l/<?) 

and constant indicator function h(6) = 1. This implies that the growth of / is not 
localized enough in the angular direction as to hope for an application of the 
Laplace method to estimate the Cauchy integrals. In other words, the second stage 
of using the saddle-point method (de Bruijn 1981, p. 77) seems to be about to 
fail. However, this is not the case here, since the whole circular contour of radius 
To = r {n) is approximately a level line of Im log z~ n f(z), not just the segment 
near the saddle point itself: in fact, from (12.6b) we get that 

a n+l/2 

Imlog/(r e ie ) -nQ = ^ sin0 + O(^ 2 ") (n ->■ 00), (12.7) 

1 — q 

which is exponentially close to zero. Note that we can arrange for r {n) £ E 
since E is built from sets of increasingly small neighborhoods of the radii of the 
zeros of /, which are located at —q~ k , k £ Nq. That is, (12.7) holds uniformly in 6. 



2 ^Note that, consistent with p = 0, we have p(r) — > as r — > 00. 
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Table 6. For f(z) = <?a( z )/ a comparison of the quasi-optimal radius 
r (w) with its asymptotic value (12.11). Note that this asymptotic value is 
not necessarily useful in practice. The value of r (n) = argmin r~ n f(r) 
was actually computed by using Matlab's f minbnd command. 



n 


A 


r*{n) 


Ko(«) 


(«/A) 2 


k(m,(«/A) 2 ) 


20 


3 


55.08575 


1.00005 


44.44444 


1.39833 


100 


15 


108.74559 


1.00000 


44.44444 


5.17900 ■ 10 11 



Hence we get 

1 r-n 1 



a « = ^-nf e- ine f(ue ie )d6 = -L- / e '( Iml °g/(^ ,e )-" e )|/(r oe ,9 )|d0 

*-7VT J — n t-TZTo J — tt 

i- /" |/(ry e )|,je.(i + o(, 2 »)), 

71~<s J— 7T 



27rrJ 

since the contribution of the odd function sin0|/(r o e ,9 )| to the integral is zero. 
Therefore, we obtain the approximation 

K«{n)=l+0{q ln ) (n -►«>), (12.8) 

whose exponentially small error term helps to understand the excellent condition 
numbers observed in Example 12.3. 

Example 12.5. We close the paper with a nontrivial example from the theory of 
random permutations. Let us denote the length of the longest increasing subse- 
quence 28 of a permutation a G S n by £(cr). The probability distribution of £(cr) 
that is induced by the uniform distribution on S n can be encoded in a family of 
exponentially generating functions <p\{z) via 

d n 
dz T > 



F(creS n :£(cr)^\)= (z) 



(A, n 6 N). (12.9) 

z=0 



Now, the seminal work of Gessel (1990) shows that <p\{z) can be expressed in 
terms of a Toeplitz determinant, 

<h(z) = det (l ]hk] (2Vz)) • (12.10) 

Since the modified Bessel functions z _fc/ ' 2 Ij : (2v / z) G No) are entire functions 
of perfectly regular growth (of order p = 1/2 and type t = 2, see Table 2), <p% 
must also be an entire function of perfectly regular growth; its order and type 
are easily inferred to be p — 1/2 and t = 2A. Likewise, we obtain that the 
Phragmen-Lindelof indicator of (p\(z) is given by 

h(8) = 2Acos(0/2) (|0|<7r). 



2 ^For instance, the longest increasing subsequence of c = (3, 7, 10, 5, 9, 6, 8, 1, 4, 2) 6 5 10 is given by 
(3,5,6,8), hence i(cr) = 4 in this case. 
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Hence, we have Cl = 1 and, since there are no zeros of <p\{z) in the vicinity of the 
real axis, to = 1 and lim n _!,oo k (w) = 1 by Theorem 10.2. This explains the very 
well behaved quasi-optimal condition numbers shown in Table 6. Theorem 8.4 
yields the following asymptotics of the quasi-optimal radius: 



However, as we can see from Table 6, this asymptotics does probably not hold 
uniformly in A and is therefore of limited practical use. Hence, one has to compute 
the value of the radius r (n) itself by numerically solving (12.2). Using these radii 
and high-precision arithmetic we were able to reproduce numerically the exact 
rational values of the distributions (12.9) for n = 15, 30, 60, 90, and 120 as tabulated 
by Odlyzko (2000), 29 who has used the combinatorial methods exposed in Odlyzko 
and Rains (2000) for his calculations. 

Remark 12.6. The numerical evaluation of <p\{z) as given by the Toeplitz determi- 
nant (12.10) turns out to suffer from severe numerical instabilities. Instead, we 
suggest to take one of the famous equivalent expressions in terms of a Fredholm 
determinant; such as the one given by Borodin and Okounkov (2000, p. 391) 



see also Basor and Widom (2000) and Bottcher (2002). Both expressions can be 
evaluated in a numerically stable way; the first using the projection method, the 
second using the quadrature method exposed in Bornemann (2010, §§5 and 6). 



I am grateful to Ken McLaughlin and Peter Miller who suggested that there 
should be a relation of the asymptotic formula (8.10) for the quasi-optimal radius 
r (n) to the saddle-point method. This has turned out to be the "missing link" 
to really understand the previously mysterious (for me at least) fact that there is 
jc^(n) « 1 for so many entire functions I had been looking to. The observation 
stated in Footnote 25 is owed to a stimulating discussion with Divakar Viswanath. 
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